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ABSTRACT 


-This  thesis  concerns  developing  a  program  that  takes  an 
aerial  photograph,  and  a  set  of  Digital  Terrain  Elevation 
Data  CDTED)  that  is  defined  over  the  area  of  the  photograph, 
and  generates  a  synthesized  view  that  represents  what  a 
camera  would  see  from  a  different  location.  The  elevation 
data  points  are  grouped  into  triangular  panels  that  are 
projected  to  the  reference  image  by  three  dimensional 
transformation  equations.  Shading  for  the  synthesized  image 
is  determined  from  the  reference  image.  The  pixels  of  the 
reference  image  that  fall  within  a  triangular  panel  are 
collected  and  averaged.  When  a  new  observer  location  is 
selected,  the  panels  are  projected  to  the  new  synthesized 
image  plane.  A  z-buffer  approach  and  a  polygon  fill 
algorithm  were  used  to  remove  hidden  surfaces  of  the 
synthesized  view. 

This  program  is  tested  on  both  artificial  and  real  data. 
Other  characteristics  and  performance  measurements  of  the 
program  are  also  analyzed  here.  The  quality  of  the  syn¬ 
thesized  image  from  real  data  was  affected  by  the  low 
resolution  of  the  terrain  elevation  data,  and  yielded  less 
desirable  results  than  could  be  expected  of  a  higher 
resolution  terrain  model. 
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I .  INTRODUCTION 


A.  COMPUTER  IMAGE  GENERATION  FROM  AERIAL  PHOTOGRAPHY 
The  main  objective  of  this  study  uas  to  develop  a 
program  that  takes  a  digital  photographic  image  and  a  file 
of  terrain  elevation  points  defined  over  that  image  as 
input,  then  produces  as  an  output  a  synthesized  perspective 
vieu.  The  synthesized  view  is  a  rotated  3-dimensional  (3D) 
perspective  representation  of  the  original  photographic 
image.  The  main  application  of  this  study  is  to  generate  a 
different  perspective  of  a  terrain  modBl .  This  may  be  used 
to  generate  different  views  that  a  pilot  of  an  aircraft 
could  expect  fallowing  different  flight  paths  through  the 
same  area.  Further  study  may  make  it  Feasible  to  generate 
synthesized  images  fast  enough  to  simulate  a  real  time  image 
display  of  a  flight  for  a  mission  briefing  or  to  be  used  as 
a  training  aid.  Another  application  could  be  For  training 
men  on  optically  guided  missiles.  With  high  resolution 
images,  a  Flight  path  through  a  battleField  could  be 
simulated  that  would  have  all  the  visual  characteristics  of 
an  actual  flight  without  the  expenditure  of  a  live  missile. 
The  generation  of  a  shaded  image  as  a  3D  picture  provides 
unique  problems  for  3D  graphic  displays.  The  data  which 
comprises  a  photographic  image  consists  of  an  array  of 
pixels,  each  of  which  has  a  defined  grey  level  cr  shade. 
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There  are  256  different  levels  of  grey  that  may  be  assigned 


to  a  pixel .  This  study  concerns  taking  a  photographic 
perspective  image  and  a  2-dimensional  C2D)  array  cf  eleva¬ 
tions  defined  m  a  grid  covering  the  area  of  the  image  as 
inputs.  T^e  grid  of  elevations,  called  the  terrain  model, 
is  geometrically  related  to  the  photographic  image  through  a 
perspective  projection  transformation  that  equates  the  world 
coordinates  of  the  elevation  points  to  the  object  coor¬ 
dinates  of  the  image.  A  synthesized  image  from  a  different 
observer  location  is  then  generated.  The  new  synthesized 
view  should  approximate  what  would  be  seen  by  a  camera  from 
the  new  observer  position. 

The  differences  between  the  original  and  synthesized 
images  will  be  affected  by  the  resolution  of  both  the 
photographic  image  and  the  terrain  models.  Higher  resolu¬ 
tion  of  the  original  models  will  result  in  a  closer  approx¬ 
imation  in  the  synthesized  image.  Another  complication  or 
ambiguity  arises  when  details  which  should  show  up  in  the 
synthesized  view  were  not  present  in  the  original  image.  A 
method  must  be  devised  so  that  it  can  fill  in  areas  which 
become  visible  in  the  synthesized  image  that  were  hidden  m 
the  original  reference  image.  The  solution  to  this  hidden 
surface  problem  is  further  addressed  in  the  discussion  cf 


the  grey  scale  referencing  a  1 gor 1 thm  .  C Ref  .  1J 
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B.  TERRAIN  ELEUAT ION  AND  PHOTOGRAPH  AS  INPUTS 

In  this  study  a  high  altitude  aerial  image  of  Moffett 
Field,  California  was  used  as  the  original  reference  photo¬ 
graph.  The  photographic  image  was  supplied  by  the  Defense 
Mapping  Agency  CDMA)  and  had  a  resolution  of  approximately  1 
meter  per  pixel .  The  terrain  model  corresponding  to  the 
reference  image  was  provided  as  Digital  Terrain  Elevation 
Data  (DTED)  by  the  DMA  and  consisted  of  elevation  points 
taken  every  second  of  a  degree  change  in  latitude  and 
longitude.  This  gives  an  approximate  resolution  of  30 
meters  per  elevation  point  in  a  north-south  direction  and  33 
meters  per  elevation  point  in  an  east-west  direction. 

The  synthesized  view  was  restricted  to  a  northerly 
direction  which  simulates  an  aircraft  flying  from  south  to 
north  with  the  image  plane  perpendicular  to  the  direction  of 
flight.  To  allow  for  different  flight  patterns  would 
require  development  of  an  algorithm  that  would  provide  for 
rotation  of  the  image  coordinates  which  is  beyond  the  scope 
of  this  study.  The  main  idea  is  to  generate  a  synthesized 
view  that  is  rotated  from  the  original  photograph  by 
approximately  90*  and  explore  the  concepts  of  the  algorithms 
required  to  do  this.  Although  speed  was  not  a  major  issue, 
the  size  of  the  terrain  model  was  limited  to  a  50  x  50  grid 
array,  cr  2500  data  points,  to  minimize  the  time  for 
synthesized  image  generation. 
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C.  ALGOR  I THn  ISSUES 

1  .  Greu  Scale  Referencing 

To  determine  the  grey  scale  values  of  the  pixels 
that  make  up  our  synthesized  view,  the  terrain  model  is 
first  divided  into  triangular  panels.  The  vertices  of  the 
triangular  panels  are  then  mapped  into  the  original  DMA 
image  using  the  perspective  projection  transformation  that 
projects  georectangular  coordinates  into  reference  image 
coordinates.  The  pixel  grey  scale  values  that  fall  within 
the  projected  triangular  panel  are  then  averaged.  The 
average  grey  scale  value  is  permanently  assigned  to  that 
particular  panel.  When  the  synthesized  image  is  constructed 
the  triangular  panel  is  mapped  to  the  new  image  view  and 
filled  with  the  assigned  average  grey  scale  value.  In  this 
way  the  sample  of  pixels  that  fall  within  the  triangular 
panel  are  mapped  from  the  original  reference  image  to  the 
synthesized  image.  This  method  of  mapping  the  triangular 
panels  to  the  synthesized  image  also  solves  the  problem  of 
assigning  a  grey  scale  value  to  hidden  surfaces  of  the 
reference  image  because  they  are  automatically  assigned  the 
value  of  the  surrounding  pixels.  Since  the  average  grey 
scale  value  for  a  triangular  panel  is  dependant  upon  the 
resolution  of  the  terrain  model,  the  grey  scale  value 
assigned  to  the  hidden  surface  will  also  be  affected 
CRef.  11. 
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The  smaller  the  triangular  panels  are,  the  smaller 
the  area  that  must  be  collected  and  averaged  in  the  original 
image.  This  means  a  much  better  synthesized  vieu  can  be 
constructed  that  contains  more  of  the  attributes  of  the 
original  image.  For  this  reason  the  resolution  of  the 
terrain  and  reference  image  model  is  very  critical  to 
obtaining  an  accurate  synthesized  vieu.  Using  the  resolu¬ 
tion  of  the  terrain  and  image  models  used  in  this  study,  the 
approximate  number  of  pixels  that  must  be  averaged  in  the 
reference  image  for  each  triangular  plane  uould  be  1/E  C30m 
x  E3m.Kl  pixel/mj  -  345  pixels.  This  is  very  coarse  and 
does  not  allou  for  optimal  generation  of  the  synthesized 
vieu  . 

S .  Hidden  Surface  Elimination 

There  are  surfaces  that  may  be  discernible  in  the 
reference  image  but  became  hidden  in  the  synthesized  vieu. 
The  z-buffer  algorithm  uas  used  to  accomplish  thB  hidden 
surface  elimination.  The  z-buffBr  is  an  array  that 
contains  the  depth  or  distance  to  the  observer  location  for 
each  pixel  that  is  to  be  visible  m  the  synthesized  image. 

As  each  triangular  plane  is  mapped  to  the  synthesized  vieu 
the  location  and  depth  of  each  pixel  uithin  the  plane  is 
determined.  The  depth  of  the  pixel  to  be  uritten  at  a 
certain  location  is  compared  to  the  depth  of  any  pixel  that 
may  have  been  previously  uritten  to  the  same  location.  If 
the  depth  or  distance  of  the  neu  pixel  to  the  observation 
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point  is  shorter  than  the  previous  pixel,  its  depth  is 
written  to  the  z-buffer  and  the  grey  scale  value  of  the 
pixel  is  placed  into  the  synthesized  image.  If  the  depth  is 
larger,  no  updating  occurs  and  a  new  pixel  is  obtained  in 
the  process.  The  z-buffer  works  very  closely  with  the  next 
algorithm  to  be  considered.  CRef.  2,  pp .  26S-267J 
3 .  Poluaon  Fill  Algorithm 

Screen  coordinates  are  generated  for  the  three 
vertices  of  each  triangular  panel  as  it  is  mapped  to  the 
synthesized  image.  Screen  coordinates  are  designated  as  IA 
and  JA  values  with  the  IA  values  representing  the  columns 
and  the  JA  values  the  raws.  The  location,  IA,JAC0,0), 
designates  the  upper  left  hand  corner  of  the  screen  and  the 
maximum  screen  coordinate,  I A , JA  C  512 , 512  3 ,  the  lower  right 
hand  corner.  An  active  edge  list  CAEL)  is  generated  by 
computing  a  line  between  each  of  the  translated  vertices. 

For  each  line,  the  IA  coordinate  corresponding  to  the 
maximum  JA  value  of  the  line,  the  amount  IA  changes  for  each 
one  unit  step  of  JA,  and  the  total  span  of  JA  are  stored 
into  the  AEL .  The  three  lines  generated  for  each  translated 
triangular  plane  will  form  another  closed  triangle.  By 
using  the  parameters  stored  in  the  AEL  the  location  of 
pixels  enclosed  by  the  translated  triangular  plane  can  be 
determined,  and  the  corresponding  array  points  within  a 
frame  buffer  are  changed  From  a  0  to  a  1. 

CRef.  2,  pp .  76-7SD 
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After  all  of  the  enclosed  pixels  have  been  marked 
within  the  frame  buffer,  the  buffer  is  scanned  row  by  row. 

If  a  value  of  1  is  found,  then  the  depth  is  calculated  for 
that  point  and  compared  with  the  depth  value  stored  m  the 
z-buffer.  If  the  depth  value  is  smaller,  that  pixel  is 
located  closer  to  the  observer  location  and  the  grey  scale 
value  for  that  pixel  is  written  to  the  synthesized  image 
file.  As  can  be  seen  the  fill  and  hidden  surface  algorithms 
work  together  to  generate  the  new  image.  The  implementation 
of  these  algorithms  are  explained  in  further  detail  later. 

In  Chapter  II  there  will  be  a  discussion  of  basic 
photographic  geometry  to  develop  an  understanding  of  the 
transformation  equations  necessary  to  map  object  coordinates 
into  image  coordinates  and  for  image  plane  rotation. 

Chapter  III  will  detail  program  considerations  based  on 
image  and  elevation  data  as  well  as  the  algorithms  used  to 
generate  the  synthesized  view.  Chapter  IU  will  discuss 
possible  ways  to  improve  the  transformation  program  and 
discussion  of  topics  for  passible  further  study.  An  outline 
of  the  program  is  contained  in  Appendix  A  that  gives  a  short 
discussion  of  each  subroutine  as  to  its  purpose,  input  and 
output,  modules  that  called,  and  modules  that  reference  the 
subroutine.  Appendix  B  contains  the  entire  3D  transforms- 


A  parallel  projection  is  a  projection  in  which  the 
projection  lines  from  the  object  to  the  image  plane  never 
converge.  When  an  object  is  viewed  by  parallel  projection, 
its  size  would  never  change  as  the  camera  is  moved  closer  or 
further  away.  In  contrast,  a  perspective  projection  has  all 
the  projection  lines  From  an  object  converge  to  a  perspec¬ 


tive  center.  A  perspective  projection  imitates  how  we  see 
things.  An  example  would  be  a  picture  of  railroad  tracks. 
The  tracks  would  appear  to  became  closer  together  when 
further  away  from  the  observation  point.  In  a  parallel 


projection  the  tracks  would  be  the  same  distance  apart  a  1  o~g 
the  entire  length.  CRef.  2,  pp .  133-1341 

Since  a  camera  is  generally  designed  to  photograph  a 
rather  large  area,  it  involves  perspective  projection.  Tie 
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camera  view  represents  ujhat  an  observer  would  see  standing 


at  the  same  location,  and  the  images  generated  are  perspec¬ 
tive  images.  This  means  that  the  equations  used  to  trans¬ 
form  the  object  space  into  the  image  space  must  be  perspec¬ 
tive  transformations. 

2 .  Image  Coordinate  Space 

The  image  plane  is  the  plane  of  the  photograph  to 
which  the  object  points  are  mapped.  It  has  a  2D  coordinate 
system  to  which  each  point  of  a  3D  object  is  translated  to 
CFig.  2.11.  The  indicated  principal  point  CIPP1  is  the 
center  of  the  image  plane  and  has  the  coordinates  of 
C x , y , 0 1 .  The  x,  y,  and  2  axis  represent  a  right  handed 
plane  and  the  perspective  center  (L)  lies  along  a  line 
parallel  to  the  z-axis;  then  a  perpendicular  line  is  drawn 
from  L  to  the  image  plane.  The  point  at  which  this  perpen¬ 
dicular  line  intersects  the  image  plane  is  called  the 
principal  point  CoJ  .  This  offset  of  the  principal  point 
from  the  IPP  is  compensated  for  in  the  transformation 
equations  by  xo  and  yo .  The  focal  length  of  the  plane  is 
defined  as  the  distance  from  the  principal  point  to  the 
perspective  center.  For  the  image  plane  coordinate  system, 
each  object  point  (A)  is  graphed  to  a  corresponding  image 
plane  point  (a)  located  at  Cxa.ya.O),  and  the  perspective 
center  or  focal  point  is  located  at  Cxo,yo,f].  CRef.  3,  pp . 
135-1363 
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In  generating  a  synthesized  view,  the  perspective 
center  may  be  placed  at  any  location  desired  with  reference 
to  the  image  plane.  It  is  therefcre  desirable  to  select  a 
point  along  the  z-axis  such  that  xo  and  yo  become  0.  This 
will  decrease  the  number  of  calculations  required  in 
generating  the  synthesized  view. 

3 .  Qbiect  Coordinate  Space 

The  observer  location  and  each  object  point  position 
is  described  m  world  coordinates,  called  georectangular 
coordinates,  of  X,  Y,  and  2.  The  center  of  the  earth  is 
given  as  C 0 , 0 , 03 ,  the  Z-axis  points  directly  to  true  north, 
the  X-axis  paints  ta  the  intersection  of  0*  Latitude  and  0* 
Longitude,  and  the  Y-axis  the  intersection  of  0*  Latitude 
and  90*  E.  Longitude.  The  principal  or  focal  point  that 
would  describe  the  observer  location  in  georectangular 
coordinates  is  XL,  YL,  and  2L .  Each  object  point  CA3 
located  in  the  object  space  is  identified  by  XA,  YA,  and  ZA 
as  shown  in  Figure  2.B.  CRef.  3,  p.  1363 

B.  inAGE  PLANE  ROTATION 

To  align  the  x,  y,  z  coordinates  of  the  image  plane  to 
the  desired  viewing  direction  requires  rotation  about  the  X, 
Y,  and  Z  axis  of  the  georectangular  coordinate  system.  In 
general  to  transform  one  30  coordinate  system  requires  a 
matrix  multiplication  of  the  form  A  “  CHIB.  The  A  repre¬ 
sents  a  vector  m  the  image  space  with  x,  y,  and  z  coor¬ 


dinates,  and  B  is  a  vector  in  the  georectangular  coordinate 


system  with  X,  Y,  and  2  components. 


This  may  be  written  as 
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where 
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This  maps  the  vector  X,  Y,  2  to  the  image  space  x,  y,  z. 

The  H  matrix  is  derived  from  the  definition  of  the  direction 
of  a  vector  ft  given  by 


CosB  j  +  Cost 


k 


(2.3) 


A  A  A 

where  i,  j,  and  k  are  unit  vectors  of  the  particular 
coordinate  system  in  which  vector  A  is  contained.  The 
quantities  c- ,  B,  t  are  the  angles  that  vector  ft  makes  with 
the  x,  y,  and  z-axis  respectively.  Since  there  are  three 
vector  components  in  the  X,  Y,  2  system,  each  one  must  make 
its  own  transformation  into  x,  y,  and  z  and  thereby  forming 


the  3x3  matrix  of  n. 


To  translate  x,  y,  z  into  X,  Y,  2 


the  inverse  of  the  M  matrix  is  taken  giving  B  “  CM]' A. 

CRef.  3,  p.  1393 

If  we  can  define  three  orthogonal  vectors  in  georect- 
angular  coordinates  as  R,  S,  and  T  that  would  describe  the 
desired  viewing  position  of  our  image  plane,  the  n  matrix  is 
easily  derived  by 


Sx 

Sy 

Sz 

Is! 

Is! 

Is! 

Rx 

Ry 

Rz 

lRl 

lRl 

lRl 

Tx 

Ty 

Tz 

HI 

in 

HI 

CE  .41 


Generally  the  image  plane  rotation  is  expressed  in  omega 
Cw),  phi  COD,  and  kappa  Ck5  .  The  uj  is  the  rotation  about 
the  X-axis,  the  9  is  rotation  about  the  Y-axis,  and  k  is 
about  the  Z-axis.  If  we  rotated  first  about  the  X-axis  by 
an  angle  of  w  radians,  the  terms  of  the  n  matrix  wculd 
equate  as  follows 

CosxX  -  CosCO* 1  -  1 
CosyY  ”  Cos(w) 

CosyZ  *  Cos(90*-w)  -  Sm(u) 

CoszY  “  CosCw  +  90*5  “  -Sin(w) 

CoszZ  *  CosCwl 

All  other  terms  equate  to  cos  90’  -  0,  therefore  the  n 
matrix  becomes 


-  A  A  ’<  A  A  A  A 


0 


0 


0  CqsCujD  Sin(ui) 


0  -Si n C  uj  1  Cos C  w  ) 


(E  .55 


Similarly  a  rotation  0  about  the  Y-axis  produces 


cos  cm  o  -sincm 


Sin  C  9  )  0 


CosC  0  D 


C  E  .  5  'J 


and  for  a  rotation  k  about  the  Z-axis 


CosC  k)  SmCk)  0 


-SinCkD  Cos  C  k  3  0  C2.7) 


By  multiplying  all  three  matrices  together  hie  derive  the 
overall  hi  transform  in  w ,  9 ,  and  k, 


Cos  C  0  )  Cos  C  k  1  CosCw)Sin(k)  ■+■  Sin  C  uj  J  Si  n  C  9  1  Cos  C  k  ) 
n  -  -Cos C 9 ) S l n C k  )  Cos  C  uj  )  CosC  k  )  -  SinCui)SinC9)SinCk} 
SinC9)  -  SinC w )CosC 9  ) 


SinCwlSinCk) 
SinC  uj  j  Cos  C  k  J 


Cos  Cui)SmC9)  Cos  C  k  ) 
Cos  C  i-j  )  S  l  n  C  9  3  S  l  n  C  k  1 
Cos  C  uj  H  Cos  C  9  3 


This  is  the  general  form  of  the  matrix  that  maps  the 
georectangular  coordinates  to  the  image  plane.  [Ref.  3,  pp 
537-6001 


The  general  form  of  the  transformation  matrix  is  used  to 
initially  map  the  elevation  terrain  model  into  the  original 
reference  image.  The  w,  9,  and  k  were  supplied  with  the 
original  DtlA  photograph  and  represents  the  rotation  of  the 
reference  image  plane  with  respect  to  the  georectangular 
coordinates  at  the  time  the  picture  was  taken.  when  we 
generate  the  synthesized  view  the  image  plane  must  be 
rotated  to  the  desired  viewing  angle  of  the  observer 
(northerly  direction  m  particular)  ,  This  means  that  a  new 
w,  0,  and  k  must  be  calculated,  or  by  defining  new  image 
plane  coordinate  axes  m  terms  of  the  georectangular 
coordinates,  we  can  calculate  the  terms  of  the  M  matrix 
directly  using  the  preceding  equations.  This  is  discussed 
further  in  the  next  chapter. 


ALGORITHM  CONSIDERATIONS 


In  this  chapter  an  in-depth  analysis  of  the  original 
image  and  terrain  data  is  presented.  This  includes  how  the 
data  was  referenced,  the  size  of  the  data  files,  how  the 
data  was  used,  and  how  the  elevation  and  image  data  compared 
with  one  another,  The  process  of  translating  from  the 
object  space  coordinates  to  image  space  coordinates  and  then 
to  screen  coordinates  is  considered,  and  the  equations  used 
are  given . 

The  referencing,  fill,  and  z-buffer  algorithms  also  are 
discussed  in  detail.  How  the  data  generated  from  these 
algorithms  is  used  and  put  together  to  produce  the  syn¬ 
thesized  view  will  be  presented.  Any  problems  that  were 
encountered  and  the  eventual  solutions  will  be  discussed  i" 
the  appropriate  section  to  which  they  pertain. 

A.  REFERENCE  I  MAGE  DATA 

The  picture  of  Naffett  Field  supplied  by  the  Defense 
Mapp:ng  Agency  CDMA),  was  4999  by  4997  pixels  in  size  and 
came  with  both  a  left  and  right  image.  Cnly  the  left  image 
was  ~sed  to  generate  the  perspective  views  m  this  study. 
Each  individual  pixel  within  the  original  image  is  desig¬ 
nated  b3  coordinates  I  and  J,  where  I  is  the  pixel  column 


and  J 


is  the  scan  row . 


The  geographic  northeast  corner  has 


the  I,  J  coordinates  of  CO,C),  and  the  southwest  corner 
C  4997 , 4999 ) . 

The  image  display  devices  used  were  capable  of  dis¬ 
playing  images  that  were  only  512  by  512  pixels  m  size, 
therefore,  the  original  image  was  divided  into  appropriate 
blacks  suitable  for  viewing  called  frames.  Each  frame 
contains  an  image  that  is  512  by  512  pixels.  The  coor¬ 
dinates  of  each  frame  has  a  four  digit  I  .Frame  and 
J_Frame  value,  and  is  further  identified  as  a  left  or  right 
CL  or  R)  image.  The  I_Frame  and  J  Frame  coordinates  are 
designated  in  multiples  of  512  which  define  the  column  and 
row  location  of  each  frame.  There  are  10  frame  columns  and 
10  frame  raws  with  assigned  coordinates  of  0000  through 
4508.  To  identify  a  particular  frame  of  the  OHA  image  one 
first  designates  whether  it  is  a  Left  or  Right  image  and 
then  give  the  I  Frame  and  J  Frame  coordinates.  As  an 
example  L05121024  would  designate  a  Left  image  From  the 
second  column  and  third  row.  The  first  frame  L00000000 
starts  in  the  southeast  corner  and  the  last  frame  L46084608 
is  in  the  northwest  corner.  The  disparity  between  the 
starting  location  of  the  frame  coordinates  and  the  I  and  J 
coordinates  of  the  original  image  must  be  compensated  for  in 
the  equations  that  are  used  to  determine  the  location  of 
individual  pixels  within  a  frame  image  as  shown  later. 
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Every  object  point  is  converted  from  its  3D  X,  V, 

and  Z  georertangular  coordinates  into  the  2D  x  and  y  image 

coordinates  using  the  A  “  ChlJB  equation  discussed  earlier. 

The  vector  from  the  perspective  center  to  each  object  point 

is  defined  by  (XA-XL),  (YA-YL),  and  (ZA-ZL).  This  vector  is 

mapped  into  the  reference  image  plane  coordinates  of  x,  y, 

« 

and  z.  Since  every  vector  in  the  image  plane  is  directed 
from  the  perspective  center  or  focal  point  to  the  image 
point  (a),  the  z  coordinate  value  is  constant  and  equal  to 
the  negative  of  the  focal  point  (-f)  of  the  camera.  Using 
these  parameters  the  equation  becomes 
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(3.1) 


uhere  K  is  a  scale  factor.  From  this  transformation  the 
following  equations  are  obtained 
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Dividing  the  last  equation  into  the  first  two  and  rear¬ 
ranging  yields 


X  "  xo  -  f 
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C3 .63 


where  x  and  y  are  the  ED  image  plane  coordinates . C Ref .  3, 
pp.  141-1423 

The  original  parameters  of  the  n  matrix  were 
calculated  from  the  ui,  9,  and  k  that  were  given  by  the  DMA 
with  the  original  photograph.  These  represent  the  physical 
position  of  the  image  plane  in  relation  to  the  georect- 
angular  coordinate  system  at  the  time  the  picture  was  taken 
The  focal  point  was  also  supplied  and  is  particular  to  the 
camera  that  was  used  to  take  the  original  photograph.  When 
the  synthesized  image  is  generated,  the  image  plane  is 
oriented  to  a  position  for  the  desired  viewing  angle,  which 
means  that  the  parameters  of  the  M  matrix  will  change  and 
must  be  recalculated.  The  steps  used  to  determine  the 
desired  orientation  of  the  image  plane  coordinates  will  be 
discussed  later  in  this  chapter. 

E .  Reference  Image  to  Screen  Coordinate  Transformation 

□nee  the  x  and  y  image  coordinates  have  been 


calculated  they  are  translated  to  I  and  J  values  of  the 


E8 


original  image.  This  is  accomplished  using  the  affine 
transform  which  represents  a  ED  into  ED  coordinate  transfor¬ 
mation.  The  equation  that  accomplishes  this  is  derived  from 


c  3 . 7 } 


where  Cl  and  CE  are  the  values  that  translate  the  image 
plane  origin  to  the  I  and  J  coordinate  system  origin.  They 
are  calculated  by  setting  I  and  J  to  0  and  solving  for  x  and 
y.  To  get  the  desired  transformation  of  the  image  coor¬ 


dinates  to  I  and  J  coordinates  the  inverse  transform  is 
taken.  CRef.3,  p.5931 
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Jm  -  kl  -1  J 


x  -  Cl 


y  -  CS 


C3  .85 


Again  the  original  J,  k,  1,  m,  Cl,  and  CE  were 
supplied  by  the  DhA  with  the  original  image.  Equation  3.8 
represents  any  general  ED  into  ED  transformation,  and  it  is 
used  both  to  translate  the  original  image  plane  coordinates 
into  the  I  and  J  coordinates  of  the  reference  image  and  to 
transform  the  image  plane  coordinates  of  the  synthesized 
view  into  screen  coordinates  of  IA,  and  JA. 

The  screen  coordinates  were  assigned  the  parameters 
IA,  which  represents  the  columns,  and  JA,  which  represents 
the  rows.  The  point  IA,JA(1,15  is  mapped  to  the  upper  left 


TiTiTi, 


corner  of  the  screen  and  I A , JAC512 , 512)  is  the  lower  right 
corner . 

The  screen  coordinates  represent  the  location  of  a 
pixel  within  a  frame.  To  convert  I  and  J  original  coor¬ 
dinates  to  I A  and  J A  screen  coordinates  one  needs  to  know 
the  particular  frame  one  is  working  in  and  compensate  for 
the  difference  in  the  starting  location  of  the  frame 
coordinates  and  the  original  image  coordinates.  This  is 
accomplished  by  using  the  following  equations, 

I A  -  Cl  -  I  Frame )  C3.S.J 

JA  -  C 4999  -  J... Frame  -  J)  C3.1G) 

The  I  .Frame  and  J.  Frame  values  must  therefore  be 


given  to  determine  the  screen  coordinates  within  a  desired 


frame.  To  allow  flexibility  in  determining  which  frame 
image  would  be  used  to  extract  the  reference  image,  an 
interactive  input  of  the  I_Frame  and  J  Frame  coordinates  was 


appropriate . 


3 .  Sunthe 


id  Image  Plane  Rotation 


For  the  synthesized  views  that  were  generated,  the 
affine  transform  parameters  Cl,  C2 ,  J,  k,  1,  and  m  that 
would  map  the  newly  rotated  image  plane  into  the  screen 
coordinate  system  were  selected.  Since  the  image  planes  m 
the  synthesized  views  were  oriented  for  an  observer  looking 
north,  the  image  plane  coordinates  were  generated  by  cal¬ 
culating  the  z-axis,  which  points  south,  from  two  georec- 
tangular  coordinate  vectors  calculated  from  two  different 
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terrain  data  points  along  the  same  longitude  line.  By 
taking  the  difference  between  the  X,  Y,  and  Z  coordinates  a 
third  vector  was  farmed  that  defined  the  image  plane  z-axis. 
The  image  plane  y-axis  was  calculated  by  using  only  one  of 
the  terrain  data  points  used  to  calculate  the  z-axis.  By 
taking  the  negative  of  the  gsorectangular  coordinates  of  the 
terrain  data  point,  a  vector  is  produced  that  points 
downward  through  the  center  of  the  earth.  This  was  the  image 
plane  y-axis.  Once  the  z  and  y  axes  are  calculated,  the 
cross  product  of  y  cross  z  was  used  to  calculate  the  x-axis. 
Figure  3.1  demonstrates  the  resulting  image  plane. 


N 


z 


Fig.  3.1  Synthesized  Image  Plane  Coordinates 

This  image  plane  coordinate  system,  from  the 
observers  perspective  at  C0,0,-f),  would  have  the  x-axis 
pointing  left,  the  y-axis  pointing  down,  and  the  z-axis 


pointing  directly  touard  the  observer.  The  screen 


coordinates  have  the  IA  axis  pointing  right  and  the  JA  axis 


pointing  down .  Therefore,  the  new  values  of  Cl  and  C2  with 
IA,  and  JA  equal  to  0  were  as  follows: 

Cl  *  Haximum  assigned  x-image  value  C3.ll) 

C2  -  0  C3 . 12) 

which  aligned  the  image  plane  x,y(o,o)  to  the  screen 
coordinates  of  IA.JACO.O),  The  J ,  k,  1,  and  m  values  of  the 
transformation  matrix  were  selected  to  scale  the  image  plane 
to  the  screen. 

Having  determined  the  three  synthesized  image  plane 
coordinate  vectors  in  terms  of  georectangular  vectors,  it  is 
relatively  easy  to  generate  the  M  matrix  parameters  using 


Sx 

Sy 

Sz 

lSl 

Is! 

Is! 

n  - 

Rx 

Ry 

Rz 

lRl 

lRl 

lRl 

Tx 

Ty 

Tz 

_  lTl 

lTl 

lTl 

were  S,  R, 

and  T 

are 

the 

C3 . 13) 


of  the  x,  y,  and  z  axes  of  the  rotated  image  plane.  This 
defines  the  new  transformation  matrix  that  will  be  usBd  to 
generate  the  synthesized  view  by  mapping  the  georectangular 
vectors  of  the  terrain  data  into  the  new  image  plane. 
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B.  TERRAIN  DATA 


The  Digital  Terrain  Elevation  Data  CDTED)  supplied  by 
the  DHA  came  as  a  rectangular  grid  of  elevation  data  points. 

Each  elevation  data  point  was  recorded  as  an  integer  value 
m  meters  above  sea  level.  If  a  particular  elevation  point 
was  unknown,  it  was  assigned  a  value  of  -33767  to  assure 
that  it  would  not  be  confused  with  any  valid  elevation  data 
points.  The  rectangular  terrain  grid  listed  elevation 
points  every  one  second  of  a  degree  change  in  latitude  and 
longitude.  The  southwest  corner  of  the  terrain  grid  was 
defined  as  being  located  at  37*  33'  47"  N.  latitude  and 
-1  EE*  05'  03"  W.  longitude.  From  this  reference  point  the 
elevation  data  was  laid  out  in  310  rows  by  339  columns.  The 
raws  represented  lines  af  constant  latitude  and  the  columns 
lines  of  .constant  longitude.  With  this  information  the 
northeast  corner  of  the  terrain  grid  was  calculated  as  being 
located  at  37*  36’  17"  N.  latitude  and  -133  01'04"  W. 
longitude . 

1  .  Data  Uenf ication 

The  first  problem  was  to  compare  how  accurately  the 
elevation  data  matched  up  with  the  original  image  data. 

This  required  taking  specific  elevation  points  that  were 
known  to  match  specific  image  points,  then  translating  the 
georectang 1 er  coordinates  of  those  elevation  points  to  the  I 
a~d  J  coordinates  for  comparison  to  the  original  image.  The 
w,  S,  and  k  for  the  n  mat:  to  transform  georectangular  to  J 


image  plane  coordinates  and  the  parameters  For  the  affine 
transform  from  image  plane  to  original  image  I  and  J 
coordinates  were  supplied  fay  the  DhA  with  the  original  image 
data  . 

To  select  specific  elevation  points  for  comparison 
required  finding  a  method  of  distinguishing  unique  elevation 
patterns  that  could  match  specific  objects  in  the  image. 

The  technique  used  was  to  visually  display  the  elevation 
data  as  an  image.  The  elevation  image  file  was  produced  by 
assigning  grey  scale  values  to  each  elevation  data  point 
with  the  lower  elevations  receiving  the  darker  shades  and 
the  higher  elevations  the  lighter  shades.  When  the  eleva¬ 
tion  image  was  displayed  as  shown  in  Figure  3.2,  a  distinct 
highway  pattern  emerged  from  which  three  intersections  could 
reasonably  be  distinguished.  The  three  intersection  points 
Cshown  as  a,  b,  and  c  in  Fig.  3.2)  correspond  to  elevated 
roads  that  crossed  one  another  in  the  original  image. 

□nee  the  elevation  points  were  selected  for  com¬ 
parison,  the  approximate  raw  and  column  of  the  elevation 
data  corresponding  to  the  center  of  the  intersections  was 
determined.  From  the  reference  point  of  the  terrain  grid 
there  is  a  1  second  change  in  latitude  and  longitudB  for 
each  row  and  column  which  allowed  the  calculation  of  the 
latitude  and  longitude  for  each  of  the  three  reference 
elevation  points.  The  latitude  and  longitude  for  each  point 
was  converted  to  georectangular  coordinates  using  a 
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conversion  program,  then  transformed  to  I  and  J  coordinates 
using  the  provided  HflA  parameters  as  explained  earlier. 


terrain  data  that  would  map  to  the  intersection,  plus  a 
small  section  of  the  surrounding  area,  was  extracted  From 
the  reference  terrain  grid.  This  smaller  set  of  terrain 
data  was  taken  from  rows  71  through  79  and  columns  172 
through  183  of  the  terrain  model  . 

To  verify  that  this  set  of  elevation  points  would 
appear  like  the  desired  reference  image,  a  line  drawing 
illustrated  in  Figure  3.3  was  created  using  a  commercial 
graphics  program  called  MOUIE.BYU  CRef.  4).  This  program 
can  generate  a  connected  line  drawing  from  a  set  of  eleva¬ 
tion  data  and  allows  rotation  as  well  as  magnification  of 
the  drawing.  To  assure  a  sharp  visual  contrast,  the 
elevation  data  was  magnified  by  a  factor  of  10  and  the 
drawing  rotated  to  a  useful  viewing  angle. 

The  results  were  not  as  desired  whic  an  early 

indication  that  the  resolution  of  the  terrain  data  may  not 
be  adequate  enough  to  generate  a  synthesized  view  that  would 
be  a  close  approximation  of  the  reference  image.  This  was 
Further  verified  when  the  synthesized  view  was  produced  at  a 
later  time.  An  artificial  set  of  image  and  elevation  data 
was  used  to  verify  that  the  transformation  program  func¬ 
tioned  as  desired  beFore  generating  synthesized  views  of  the 
original  reference  image.  The  artificial  elevation  points 
mapped  into  the  artificial  image  exactly  way  as  the  original 
elevation  data  would  map  to  the  original  reference  image. 

The  artificial  reference  image,  shown  in  Figure  3.4, 
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uas  of  a  small  square  structure  resting  on  top  of  a  much 
larger  square  structure.  Selecting  a  large  object  for  an 
artificial  reference  image  would  minimize  the  effects  of  the 
low  resolution  of  the  elevation  data.  This  produced  more 
reasonable  synthesized  images  that  demonstrated  the  perspec¬ 
tive  transformation  more  accurately.  The  results  of  the 
transformation  will  be  discussed  further  on  in  this  chapter 
after  considering  some  of  the  algorithms  used  to  generate 
the  synthesized  view. 

C.  SPECIFIC  ALGOR  I  THUS 

1 •  Image  Referencing  Algorithm 

Due  to  the  resolution  mismatch  between  the  reference 
image  and  the  terrain  data,  a  smaller  subset  of  the  terrain 
model  uas  selected.  This  uould  allow  the  transformation  of 
smaller  and  more  distinct  images  that  could  show  the  effects 
of  the  program  better.  The  desired  terrain  elevation  points 
are  extracted  from  the  larger  reference  terrain  grid  by 
interactively  selecting  the  proper  rows  and  columns  that 
define  those  particular  points.  The  program  allows  up  to  a 
50  by  SO  array  of  elevation  points  to  be  extracted.  This 
size  uas  chosen  to  limit  the  time  required  to  generate  a 
synthesized  view.  From  the  rows  and  columns,  the  latitude 
and  longitude  is  tabulated  For  each  point,  and  is  then 
converted  to  georectangular  coordinates.  The  georectangular 
coordinates  are  written  to  a  file  in  row  order  from  left  to 


right  and  from  bottom  to  top.  The  first  set  of  coordinates 
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match  the  scuthujest  elevation  point,  and  the  last  set  of 
coordinates  the  northeast  elevation  point. 

The  georectangular  coordinates  of  each  terrain  data 
point  is  translated  to  I,  and  J  original  image  coordinates 
using  the  camera  parameters  supplied  by  the  DMA.  They  are 
further  processed  to  IA  and  JA  screen  coordinates  using  the 
I  Frame  and  J  Frame  of  the  desired  reference  image.  The  IA 
and  JA  coordinates  derived  from  the  elevation  points  uiil 
now  correspond  to  the  IA  and  JA  coordinates  of  the  reference 
image.  Any  four  adjacent  elevation  points  will  define  a 
rectangle  which  is  divided  into  two  triangular  panels  by 
defining  a  line  connecting  two  opposite  corners.  Once  the 
triangular  panels  of  the  elevation  points  are  defined  m 
terms  of  IA  and  JA  coordinates,  the  corresponding  pixel  grey 
scale  values  of  the  reference  image  that  fall  within  the 
same  screen  coordinates  of  the  triangular  panel  are  col¬ 
lected  and  averaged.  As  seen  in  the  example  of  Figure  3.5, 
the  calculated  average  grey  scale  value  is  placed  into  a 
file  along  with  the  three  specific  elevation  points  that 
make  up  the  triangular  panel.  This  procedure  is  repeated 
until  the  entire  terrain  grid  has  been  processed.  Each 
triangular  panel  represents  a  sample  or  average  grey  scale 
value  of  the  reference  image. 

When  the  latitude,  longitude,  and  elevation  cf  the 
new  observation  point  is  input  into  the  program,  a  new  XL, 


YL,  and  ZL  is  calculated  that  corresponds  to  the  new 


perspective  center.  As  explained  earlier,  the  image  plane 
is  rotated  to  the  desired  viewing  angle,  which  changes  the  n 
matrix  parameters  as  well. 
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Fig.  3.5  Triangular  Panel  Data  File  Structure 

With  these  new  parameters  For  the  perspective  transformation 
equations,  the  georectangular  coordinates  of  the  terrain 
grid  are  once  again  run  through  the  perspective  transfor¬ 
mation,  and  the  new  IA  and  JA  coordinates  are  calculated  for 
each  point.  These  new  IA  and  JA  coordinates  represent  the 
transformed  elevation  points  as  seen  from  the  new  observer 
location.  The  next  step  is  to  take  the  same  three  elevation 
points  that  farmed  a  triangular  panel  in  the  original 
transformation,  map  them  into  the  synthesized  image  file 


using  the  new  IA  and  JA  coordinates,  and  then  fill  them  m 


In  the  perspective  transformation  of  the  terrain 
grid  into  the  synthesized  image  plane,  seme  of  the  tri¬ 
angular  panels  become  partially  or  entirely  hidden  by  other 
panels.  To  determine  which  pixels  will  be  visible  and 
therefore  written  to  the  synthesized  image  and  which  pixels 
are  hidden,  the  z-buffer  algorithm  was  used. 

When  the  new  observer  location  is  input  into  the 
pregram,  the  XL,  YL,  and  2L  georectangular  coordinates  are 
tabulated.  Using  the  georectangular  coordinates  previously 
calculated  far  each  of  the  terrain  grid  elevation  points, 
the  distance  or  depth  from  the  observer  location  to  the 
elevation  data  paints  can  be  calculated  by  the  following 
equation 

Depth  -  Square  rootC CXL-X)*  ♦  CYL-Y)*  +  (ZL-Z1 1 C3.14) 
The  depth  of  the  pixels  within  a  triangular  panel 
were  calculated  by  using  the  normalized  plane  equation  that 
defines  the  plane  of  the  triangular  panel.  The  normalized 
plane  equation  in  3D  space  is  given  by 
ax  +  by  ♦  cz  •  -1  C  3 . 1 5  j 

To  solve  for  the  coefficients  a,  b,  and  c,  the  three  eleva¬ 
tion  points  that  specify  a  triangular  panel  are  used,  and 
the  <,  y,  and  z  are  replaced  with  IA,  JA,  and  the  DeptH  of 
each  elevation  point.  If  the  three  points  are  (IAl,  JA1 , 


Depthl) ,  (IA2,  JA2,  Depth2),  and  CIAS,  JA3,  Depth3),  then  in 
matrix  farm  we  have  the  fallowing 


IA1  JA1  Depthl 
IA2  JA2  Depth2 
IA3  JA3  Depth3 

Salving  far  the  coefficients 
2083 


a,  b,  and  c  we  have  CRef.  2,  p. 


I A 1 

JA1 

Depthl 

-1 

-1 

IA2 

JA2 

Depth2 

-1 

I  A3 

JA3 

Depth3 

-1 

The  inverse  af  the  elevation  point  matrix  is 
determined  by  calculating  the  adjoint,  which  is  the  trans¬ 
pose  of  the  cofactar  matrix,  and  multiplying  each  term  by 
the  reciprocal  af  the  determinant  CRef.  5,  p.  A-153  .  The 
depth  of  each  pixel  within  a  tranformed  triangular  panel  can 
now  be  determined  by  using  the  IA  and  JA  of  the  pixel  in  the 
following  equation. 

Depth  -  -  Cl  +  a  CIA)  +  b  C JA 3 )  /c  C3.1B) 

Each  triangular  panel  will  have  different  a,  b,  and  c 
coefficients,  therefore  the  pixels  within  each  triangular 
panel  will  have  a  different  depth  than  those  from  another 
panel  . 

The  z-buffer  is  a  two  dimensional  array  that  is  the 
same  size  as  the  synthesized  image  of  512  by  512  pixels. 
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LJhen  the  IA  and  JA  coordinates  of  a  pixel  is  determined,  the 
depth  is  calculated  and  then  compared  to  any  previously 
written  depth  in  the  z-buffer  at  that  same  coordinate 
location.  If  the  depth  is  smaller,  it  means  that  the  pixel 
is  closer  to  the  observer  and  would  cover  the  previously 
written  pixel.  The  depth  of  the  new  pixel  will  then  replace 
that  of  the  previous  pixel.  The  assigned  grey  scale  value, 
which  is  determined  from  the  projected  triangular  panel  in 
which  the  new  pixel  is  contained,  is  written  to  the  syn¬ 
thesized  image  at  the  calculated  IA  and  JA  coordinate 
location.  If  the  new  pixel  depth  is  larger  then  the 
previous  pixel  depth,  that  means  the  new  pixel  is  farther 
away  from  the  observer  and  would  be  covered  by  the  previous 
pixel.  In  this  case  no  action  is  taken  and  the  next  pixel 
is  processed.  This  procedure  is  continued  until  all  the 
pixels  within  each  translated  triangular  panel  has  been 
processed.  CRef.  E,  pp .  E65-E67J 

The  IA  and  JA  coordinates  of  the  pixels  in  the 
synthesized  image  are  calculated  from  each  transformed 
triangular  panel  using  an  active  edge  list  CAEL),  J  Bucket, 
and  frame  buffer.  The  IA  and  JA  screen  coordinates  of  the 
three  elevation  points  that  define  a  triangular  panel  are 
used  to  generate  the  parameters  of  the  AEL .  Three  lines  are 
f armed  that  connect  each  of  the  three  translated  elevation 
points  and  are  identified  as  line  1,  S,  and  3.  A  line  is 
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defined  between  any  two  points  from  which  you  can  determine 
the  I A  associated  with  the  maximum  JA  of  the  two  points  by 
I A  I NCPT  -  C  IA  value  of  the  maximum  JA  C3.19) 

How  much  IA  changes  for  a  one  step  change  in  JA  is  given  by 


DELTA  I A 


IACPomt  2?  -  IACPomt  1 } 


JACPoint  2)  -  JACPamt  1) 


C3 .20) 


which  is  the  inverse  slope  of  the  line  between  the  two 
points.  The  span  of  JA  between  two  points  is  given  by 

DELTA, JA  -  JAC  PointE )  -  JACPoint  1)  C3.21) 

For  each  line,  the  IA  coordinate  that  corresponds  to 
the  maximum  JA  value  of  the  line,  the  amount  IA  changes  for 
each  one  unit  step  of  JA,  and  the  total  span  of  JA  is 
determined  and  stared  into  the  AEL .  After  the  line  para¬ 
meters  are  placed  into  the  AEL,  the  line  identification 
number  is  put  into  the  J  Bucket ,  which  is  a  one  dimensional 
array  of  size  512,  at  the  same  location  as  the  maximum  JA  of 
the  line.  In  this  way  the  JBucket  acts  as  a  pointer  to  the 
lines  stored  in  the  AEL.  If  a  previous  line  number  has 
already  been  written  to  that  location,  then  the  line  ident¬ 
ification  number  of  the  new  line  is  placed  into  the  AEL  and 
is  linked  to  that  line  already  residing  in  the  J, Bucket.  An 
example  is  shown  in  Figure  3.6  where  line  #1  is  referenced 
by  the  J  Bucket  and  line  #2  is  linked  to  line  #1.  The  IA 
and  JA  coordinates  of  the  pixel  locations  to  be  mapped  to 


the  synthesized  image  are  contained  within  the  three  lines 


whose  parameters  are  contained  m  the  PEL.  If  a  line  is 
horizontal,  the  IP  and  JP  coordinates  of  that  line  are 
located  between  the  other  two  lines  and  therefore  does  net 
need  to  be  written  to  the  PEL.  [Ref.  2,  pp . 75-703 

J  Bucket  PEL  For  Line  1 


PEL  For  Line  2 

Fig.  3.6  Pctive  Edge  List  Storage 


The  J _Bucket  is  scanned  from  its  maximum  to  minimum 
coordinate  value.  If  the  J  Bucket  contains  a  line  pointer, 
the  parameters  for  that  line  are  retrieved  from  the  PEL  as 
well  as  those  of  any  line  it  is  linked  to.  Since  we  have 
defined  a  triangle  there  will  always  be  two  lines  to  work 
with.  From  the  parameters  of  the  two  lines,  the  IP  coor¬ 
dinates  that  fall  between  them  is  calculated  for  each  JP 
scan  line.  P  scan  line  is  all  the  columns  CIP  coordinates) 
along  a  particular  row  (JP  coordinate).  Ps  the  JP  coor¬ 
dinate  is  decremented  by  one,  from  maximum  to  minimum,  the 
J  Bucket  is  checked  to  see  if  a  new  line  has  been  added,  and 
then  the  corresponding  IP  for  each  line  is  determined  for 
that  JP  value.  Those  IP  values  and  any  that  fall  between 
them  are  mapped  to  the  frame  buffer  . 


The  frame  buffer  is  a  515  by  515  array  that  is 
initialized  to  0.  As  the  IA  and  JA  values  are  mapped  into 
the  frame  buffer  the  coordinate  locations  matching  the  IA 
and  JA  values  are  changed  from  0  to  1 .  This  process 
continues,  and  each  time  the  JA  scan  line  is  decremented, 
the  DELTA  JA  parameters  for  each  of  the  two  lines  are  also 
decremented  and  the  J  Bucket  is  checked.  If  the  J  Bucket 
contains  another  line  pointer,  then  the  line  it  references 
will  replace  the  line  whose  DELTA  JA  has  decreased  to  0. 

When  finished,  the  frame  buffer  will  have  recorded  the  IA 
and  JA  coordinates  for  every  pixel  location  within  the 
translated  triangular  panel.  The  frame  buffer  is  then 
scanned,  and  if  a  particular  IA  and  JA  coordinate  location 
contains  a  value  of  1,  then  the  depth  for  a  pixel  located  at 
those  same  coordinates  is  calculated  and  compared  to 
previously  tabulated  depths  in  the  z-buffer  as  explained 
ear  1 ler . 

This  process  is  known  as  a  polygon  fill  routine  that 
utilizes  the  z-buffer  algorithm  for  hidden  surface  elimina¬ 
tion.  If  any  part  of  a  triangular  panel,  after  going 
through  the  perspective  transformation,  should  map  outside 

I 

j  the  synthesized  image  coordinate  boundary,  the  entire  panel 

is  discarded  and  the  next  panel  is  processed.  This  is  not  a 
satisfactory  solution  and  could  have  been  corrected  by 
implementing  a  clipping  routine  that  would  allow  partial 
triangular  panels  to  be  mapped  to  the  synthesized  image. 

i 
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However,  due  to  time  constraints,  a  clipping  routine  was  net 
incorporated  into  the  3D  transformation  program. 


□.  SUHMARY 

Usi^g  the  artificial  reference  image  from  Figure  3.4  and 
an  artificial  terrain  grid  that  mapped  to  the  same  image, 
two  synthesized  views  were  generated.  The  first  synthesized 
view  shown  in  Figure  3.7  was  from  an  observation  point 
located  at  37*  22'  30"  N.  latitude,  -122*  01'  59"  W. 
longitude,  and  an  elevation  of  110  meters.  Figure  3.8 
exhibits  the  next  synthesized  view  that  was  produced  from  an 
observation  point  of  37*  20'  0”  N.  latitude  while  main¬ 
taining  the  same  longitude  and  observer  height  as  before. 
This  simulates  viewing  the  object  from  further  away.  As 
expected  of  a  perspective  view  the  object  appears  smaller. 
The  near  view  also  demonstrates  the  perspective  relationship 
by  the  apparent  taper  from  front  to  back.  A  third  perspec¬ 
tive  view  is  depicted  in  Figure  3.9  from  close  in  and  at  a 
higher  elevation.  The  observer  location  was  from  37*  23’ 

30  1  '  N.  Latitude,  -122*  01'  59''  W.  Longitude,  and  a  height 
of  150  meters.  This  demonstrates  the  visual  effect  of  not 
displaying  partial  triangular  panels  in  the  image  and  why  a 
clipping  routine  is  needed. 

The  actual  reference  image  is  shown  in  Figure  3.10  from 
which  the  synthesized  view  displayed  in  Figure  3.11  was 
generated.  In  comparing  the  reference  image  to  the  syn¬ 
thesized  view  there  is  little  resemblance.  A  closer 
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synthesized  view  is  seen  in  Figure  3. IS  which  gives  a 
partial  representation  of  the  reference  image.  The  lack  of 
detail  in  doth  shading  and  the  depiction  of  features  can  be 
attributed  to  the  poor  resolution  of  the  terrain  data  (ap¬ 
proximately  55  meter  resolution)  as  compared  to  that  of  the 
reference  image  data  Cl  meter  resolution). 

To  improve  the  quality  of  the  synthesized  view  the 
terrain  data  must  be  of  a  higher  resolution.  Having  more 
elevation  data  points  over  a  given  area,  the  fewer  number  of 
reference  image  pixels  one  must  collect  and  average  for  a 
triangular  panel  .  The  synthesized  image  will  then  maintain 
a  closer  approximation  to  the  various  shades  of  the  refer¬ 
ence  image.  The  physical  shape  of  an  object  also  suffers 
from  poor  terrain  data  resolution.  The  perspective  trans¬ 
formation  forms  straight  lines  between  translated  elevation 
points.  This  causes  object  distortion  if  the  elevation 
points  do  not  fall  exactly  along  the  boundaries  of  the 
object.  Os  an  example,  if  a  square  box  in  the  reference 
image  had  only  one  elevation  point  defined  on  its  surface  at 
the  center  of  the  box,  then  the  synthesized  image  can  not 
reproduce  the  corners  and  edges  of  that  box,  and  it  would 
appear  distorted.  For  these  reasons  the  closer  the  resolu¬ 
tion  of  the  terrain  data  to  the  resolution  of  the  reference 
image,  the  closer  the  synthesized  view  resembles  the  refer¬ 
ence  image  in  shading  and  shape. 


Fig.  3.9  Artificial  Synthesized  Image  3 


Fig.  3.10  Reference  Image 


Fig  3.11  Synthesized  Image  1 


The  30  transformation  program  was  developed  to  allow 
tracing  the  flow  of  data  easily.  riuch  of  the  data  was 
written  to  files  so  that  it  could  be  printed  out  and 
studied.  This  method  of  data  storage  required  certain  fil 
to  be  read  many  times  as  the  synthesized  image  was  gener¬ 
ated.  The  program  execution  would  have  been  faster  if  the 
data  had  been  stored  in  arrays  that  are  easily  passed 
between  the  various  modules.  Another  consideration  that 
would  increase  speed  would  be  to  decrease  or  eliminate  the 
interactive  input  required.  This  could  be  accomplished  by 
extracting  the  desired  terrain  data  into  a  file  and  sepa¬ 
rating  the  associated  reference  image  before  program 
execution.  This  would  require  longer  set  up  time  but  would 
decrease  the  amount  of  data  manipulation  required  by  the 
program  and  increase  its  overall  speed. 


IU .  CONCLUSIONS 


A.  GENERAL 

The  3D  computer  image  transformation  from  a  photographic 
image  mas  a  difficult  task  to  achieve  satisfactorily.  The 
main  objective  was  to  develop  a  program  that  takes  a  grey 
scale  photographic  image,  a  set  of  elevation  data  points 
defined  over  that  image,  and  generates  a  rotated  synthesized 
perspective  view,  This  goal  was  realized,  but  some  areas 
still  need  improvement.  The  quality  of  the  synthesized  view 
is  Judged  on  how  well  the  grey  scale  values  of  pixels  match 
those  of  the  original  image  and  how  closely  the  translated 
objects  resemble  the  desired  structure. 

The  resolution  of  the  terrain  model  as  compared  to  the 
resolution  of  the  reference  image  data,  extensively  affects 
the  synthesized  image  quality.  Depending  on  the  contents  of 
the  reference  image,  the  required  resolution  of  the  terrain 
model  mill  vary.  If  the  image  is  of  open  country,  the 
distance  between  elevation  points  may  be  large  and  the 
result  may  not  suffer  unacceptable  degradation  in  the  syn¬ 
thesized  view.  If  the  image  is  of  a  city  that  has  many 
small  distinct  objects  such  as  buildings,  then  the  resolu¬ 
tion  of  the  elevation  points  must  be  higher.  Given  a  set  of 
elevation  data  points,  the  question  to  be  resolved  is  hem  to 
make  the  synthesized  pixels  relate  to  the  reference  image 
better,  and  thus  improve  the  quality  cf  the  synthesized 
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views7  This  question  and  alternative  ways  to  improve  speed 
and  program  flexibility  will  be  discussed. 

B.  GREY  SCALE  CORRELATION 

There  are  a  few  possible  methods  to  improve  the  shading 
of  the  synthesized  view  to  better  match  that  of  the  original 
image.  The  translated  triangular  panels  form  very  distinc¬ 
tive  lines  or  boundaries  between  areas  of  different  shades. 
It  may  be  desirable  to  blend  these  boundaries  by  sampling 
the  pixels  along  both  sides,  replacing  those  pixels  with  an 
averaged  value.  Another  possibility  would  be  to  implement  a 
Gouraud  or  Phong  shading  algoritm  CRef,  3,  pp .  323-3301. 

This  would  help  smooth  the  shade  transition  across  the 
boundary  and  result  in  a  more  smooth  appearance. 

Another  method  to  improve  grey  scale  correlation  would 
be  to  develop  a  way  to  divide  the  triangular  panels  into 
smaller  triangular  areas  before  referencing  the  original 
image.  By  having  smaller  triangular  panels,  smaller  areas 
of  the  reference  image  are  sampled  resulting  in  a  better 
approximation  in  the  synthesized  view. 

Only  the  left  image  of  a  stereo  photographic  pair  of 
images  was  used  in  this  study.  The  right  image  may  contain 
grey  scale  information  of  surfaces  hidden  m  the  left  image 
view  or  vice  versa.  By  sampling  both  left  and  right  images 
and  complementing  them,  some  ambiguities  may  be  resolved. 
These  suggestions  will  increase  the  number  of  calculations 
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required  to  generate  the  synthesized  view  but  may  improve 
the  quality  of  the  synthesized  view. 

C.  PROGRAM  SPEED  AND  FLEXIBILITY 

Although  speed  was  net  a  prime  consideration  in  this 
study,  it  would  be  desirable  to  generate  synthesized  views 
as  quickly  as  possible.  To  determine  which  subroutines 
consumed  the  most  GPL!  time  the  program  a.as  mentored  while 
running  the  artificial  reference  data  set.  Table  1  shows 
the  results  of  the  CPU  time  used  and  the  percentage  of  the 
total  time  far  each  subroutine  called  by  the  main  program. 

If  a  subroutine  calls  another  subroutine  that  time  is 
included  m  the  calling  routines  CPU  time.  The  subroutines 
that  required  interactive  input  were  not  measured.  There¬ 
fore,  the  total  CPU  time  cannot  be  measured  precisely. 
However,  the  results  obtained  represent  reasonable  estimates 
and  demonstrate  where  the  focus  should  be  on  improving 
overall  speed. 

Some  suggestions  have  already  been  discussed  on  how  to 
improve  the  speed  of  the  program,  but  other  methods  also 
exist.  Preconditioning  the  input  data  to  eleminate  interac¬ 
tive  input  and  using  arrays  instead  cf  files  were  two 
methods  already  presented.  The  z-buffer  is  accessed  several 
times  during  synthesized  view  generation.  If  the  z-buffer 
could  be  implemented  m  hardware,  the  time  required  for 
processing  of  polygon  fill  and  hidden  surface  elimination 
routines  would  be  improved. 
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TABLE  1 

CPU 

TINE  CONSUnPTION 

C IN  SECONDS) 

SUBROUTINE 

CPU  TINE 

PERCENTAGE  C  % ) 

TER  CROP 

14.11 

9.051 

READ  I  flAGE 

3.23 

2.072 

TER  INTRP 

0.04 

0.026 

REAL  EL 

0.07 

0.045 

TER  DUS 

0 .69 

0.443 

in  REF I J 

0.39 

0.250 

in  REFAUG 

1 .10 

0.706 

AFFIN 

0.04 

0.026 

NEW  IJ 

0.90 

0.513 

NODE  DEPTH 

0.63 

0.404 

FILL 

134 . BO 

96.467 

TOTAL 

155.9 

100.0 

The  frame  buffer  is  used  and  accessed  in  two  different 
areas  of  the  program  for  each  triangular  panel  translated. 

It  may  be  possible  to  eliminate  this  buffer  by  processing 
each  pixel  as  its  IA  and  JA  coordinates  are  determined, 
instead  of  storing  that  information  into  the  frame  buffer 
for  later  processing.  Other  polygon  fill  routines  such  as 
seed  fill  algorithms  or  using  fence  registers  may  be  faster 
than  the  edge  fill  routine  used  in  this  program  CRef.  2,  pp . 
00-063  . 

For  program  flexibility  it  would  be  desirable  to  be  able 


to  adjust  the  image  plane  of  the  synthesized  view  to  any 
desired  angle.  The  program  limits  the  image  plane  to  a 
northern  direction.  To  make  the  synthesized  image  plane 


adjustable  would  require  developing  a  method  for  rotating 
the  image  plane  coordinates  m  terms  of  the  georectangular 


coordinates,  This  would  improve  the  3D  transformation 
program  and  extend  its  usefulness.  Another  program  improve¬ 
ment  would  be  the  incorporation  of  a  clipping  routine  to 
improve  the  appearance  of  partial  synthesized  views  as 
discussed  m  previous  sections. 

The  ideas  used  to  develop  this  program  were  contrived 
from  fundamental  concepts.  flany  areas  were  discussed  that 
could  improve  or  enhance  the  basic  implementation  of  the 
program.  The  results  that  were  obtained  are  encouraging  and 
could  easily  be  used  as  the  basis  for  further  study. 


^a»»ww' 


A 

■ 

I 

’S 

\ 

> 

» 

"i 

► 


„ nrKnnrr. nr  wvv u*VY\-vvY\rr W  V> V» >-*  ST*  V  -j.  -j-  v  -r  .-w  *--  r 


CfSjlf  available  id  DTlC  Hoes  not 
g*nnii  fullj  legibly  reproduction 


PPPEME-X  ^ 

ppgcp^’-t  ci  iM^cipy 


1  .  Pr22P^K‘  3  3 

3.  Fw^^ti'c^s  csr~Fwr'fr'2d 

Takes  =n  inage  ard  slsvatic"  data  file  a~d  e  tree: 
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The  file  containing  the  reference 
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An  array  containing  the  extracted 
elevation  points. 


d.  „a**.ng  routines 

Tppf.j  3  3 , 

e.  Tailed  routines. 
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a-  d  Zou-ters  . 
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e  . 


X  IMA  : 

The  image  plane  x-axis 

coord i ~  a  t  e 

value  of  the  extracted 

data  points . 

elevat to" 

y  t  yft  . 

The  image  plane  -.-axis 

coordinate 

value  of  the  extracted 

points  . 

elevatic~ 

Zal-irg  rcuti 

I  -'  PEFIJ. 

nas 

Tailed  routi° 

None  . 

es 

Routine  parameters 

hil,  Ml? ,  r i 3 

• 

MSI ,  M22,  M23 

M31,  M32 ,  ~33 

:  The  parameters  of  the  M 

matri  . 

T-|  *  |HM 

The  denominator  of  the 

tion  equations. 

transf or^a 

‘CUT  I  IE  XY2IJ 

F.ncticns  per 

formed 

Eo"ver ts  the 

image  plane  coordinates  of 

tl  ^  2  S  .  :  ~ 

ial  i~age 


:ted  elevation  data  paints  ta  the  I  and  J  art 
■  d  .-ates . 

Input 

XIM(5;  The  .mage  plane  .-a  is  coordinate 

values  of  the  extracted  elevation 
data  paints. 


I  ft  and  Jft: 

The  arrays  that  pent 

fere^cs  image  seres’- 

the  e,:traoted  elevat 

I^3£: 

pCi--S . 

The  array  that  ccnta 

reference  image  grey 

I  END!"!: 

The  number  of  rous 

elevation  data. 

I  t  iu  J  : 

The  number  of  oolum^ 

extracted  elevation  data. 

□utput 

NDCE . D^T :  The  Tile  that  contains  the  eisva 

t ion  points  that  make  up  a  tri¬ 
angular  panel  and  its  associated 
reference  grey  scale  value. 

Tailing  routines 

TFo-j  3  3  . 

Tailed  routines 

hone  . 

Routine  parameters 

SLCPE :  The  slope  of  the  diagonal  l;ne  t 

separates  twe  rectangle,  defined 
four  elevation  data  points,  ."to 
tuo  tna"3ilar  panels. 
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13.  SUBPOUT  I ME  EXT  OPIET 


F.":t::rs  performed 


Determines  the  M  matrix  parameters  used  m  the 


itation  cf  the  mage  plane  far  the  synthesized  vieu.. 


:rogram  defines  the  three  mage  plane  coordinates  m  ter 


;ecreatangui ar  coordinates. 


t.  Input 


IENEN: 


The  number  of  raus  in  the  extra 


elevation  data  points. 


I  ENDN 


The  number  of  columns  in  the 


extracted  elevation  data  points 


Output 


ril,  TIE,  HI3:  First  row  coefficients  of  the 


transformation  matrix. 


031,  ri33 ,  T23:  Second  rou  coefficients  of  the 


ran.sf ormat ion  matrix. 


T31  ,  "'33,  ri  3  3  :  Third  rou  coefficients  of  the 


transformation  matri: 


lailing  routines 


Me'u  1  . 


Sailed  routines 


Ujt;ne  oarameter 


nAGN 


The  magnitude  c f  the  synthesized 


is j  i mage 

p  lane 

: : -a. 

:  i  s  . 

PAGN  V: 

The  nagn 1 1_ 

de  cf 

the 

Synthesized 

vise  it, age 

plane 

y-a: 

:  1  s  . 

» /S  r-~  f  ~7 

^  0-  '  «  l—  - 

The  -t a g n 1 1 _ 

da  c f 

*-  V- 

synthesized 

vieu  mage 

plane 

Z-ax: 

: is  . 

X,  V,  and  2: 

Georectangu 

1  a  r  co 

or  d  i 

nates  cf  i'nE 

extracted  elevation  data  points. 

X  CORD:  Array  that  stores  the  gecrec- 

t angular  X  coordinates  of  the 
extracted  elevation  data  points. 

V  CCPD :  Array  that  stores  the  gecrec- 

tangular  V  coordinates  of  the 
extracted  elevation  points. 

2  OCPD:  Array  that  stcrss  the  gecrec- 

tang^iar  2  coordinates  of  the 
extracted  elevation  data  pcmts. 

'  -JECX,  X  CEO  1 , 

and  X  PEC2:  The  Y  .  ard  2  gecracta^g-lar 

laordi'-ates  of  the  synthesized 
i nape  plana  /.-axis, 

i  111.!,  t’  v E 0  , 

and  VEC2 :  The  X,  t  ,  anp  z  gecrecoang_lar 

coordinates  of  the  Sy -;ues *  zed 
inage  plane  y-a.us. 


Z  UECX,  Z  UECY, 


and  Z  O'ECZ :  The  X,  Y,  a^d  Z  gecrscta'-g-lar 

oordmates  of  the  s:j'-tNes.zed 
I’s-js  plane  z-ax:s. 

ITCT:  The  total  number  of  extracted 

elevation  points. 

IP  a^d  IP:  Counters. 


□  JUNE  OFF  I N 
unctions  performed 

Assigns  or  caloulates  the  ooeffioients  to  be 
,  m  the  affine  transform  from  image  coordinates 
oordinates  of  the  synthesized  visu. 

I ’-put 
None  . 

p  L-.  ^ 


Ol , 

a?  , 

si, 

E2 

and 

T-  i 

L) 

C2: 

The  affi’-e  transform  ooeff.o 

for  the  synthesized  vieu  image 
plans  ooord i-ates  to  scrs2r 
coordinates . 

Calling  routines 
TPPN  3  C . 


X 


VI,  and  21: 


he  gscrectangu  1  ar  coordinates  of 


neui  observer  location. 

FOCUS:  Assigned  fecal  length  for  the 

synthesized  vieu. 


Cal. mg  routines 
CRAM  3  0  . 

Called  routines 
CMS2XV2 . 

Routine  parameters 


f  icne 


IE.  SUBROUTINE  NEU  IJ 

a.  Functions  performed 

Calculates  the  new  I  A,  and  JA  screen  coordinates  of 
the  extracted  elevation  data  points  using  the  new  observer 
location  data . 

b .  I nput 

XI,  VI  and  21:  Ceorectangul3r  coordinates  of  t'-e 

new  observer  location. 

FOCUS :  Ass^g^ed  focal  length  for  the 
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IENDN:  The  number  cf  cci^nns  : 

extracted  e I evat: :r  jata  pei-t: 

Output 

I  A:  Tus  arrc3  aa^tsirirg  the  Slavs' 

data  I  A  scree1-  cccrdi-ates  cf 
synthesized  vis.  , 

_T-:  The  array  cc-tar-i-s  the  alEva1 

data  J~  screen  cccrdi^atss  cf  1 
synthesized  '.tsi  . 

Calling  rectmes 

TP»N  3  D  . 

Called  routines 

y y 3  t  j  EXT  CP I EN 

Rautine  pararrsters 

f "  1 1 ,  ~12,  and 

ri3:  First  rcu  ccefficients  cf  tHe 

trarsf mat ion  matrix. 

-‘21,  ”22,  and 
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YIMA  : 


Image  plane  y  coordinate. 


□ENOn : 


I  TCP 


IF  : 


Denominator  of  the  transformat i : 
matrix  . 

Total  number  of  extracted  eieval 
data  points. 

Counter . 


SL'ERCC'T  I NE  NODE  DPTH 

a.  Functions  performed 

Calculates  the  distance  from  each  transformed 
vation  data  point  to  the  new  observer  location. 

b .  I nput 

XI,  Yl,  and  21:  The  georectangu 1 ar  coordinates  of 

the  new  observer  location. 

lENCh:  The  number  of  rows  m  the  extracts 

elevation  data  points. 

I  E\D*J :  The  number  of  columns  m  the 

extracted  elevation  data  points. 


DEPTH 


Arr  aa  of  the  distances  from  the 
transformed  elevation  data  points 
to  the  "eu  observer  lccat*cn. 


c.  balling  r^utres 
TRAN  3  D. 

e.  Cs.led  routines 
None  . 

f.  Routine  parameters 


h  CuEF ,  B  CGEF 


and  C  CDEF: 


'hs  csef f icierts  of  the  plan 


equation . 


I  GREY : 


Grey  scale  value  of  a  tr;arg 


►  j n r> r  £5  M  E r  c 


panel  . 


Turner ical  dssiqnat.c^  of  the 


elevation  data  points  that  n 


triangular  panel. 


Minimum  13  screen  coordinate 


the  three  elevation  points, 


Mar. imum  screen  coordinate 


three  elevation  points. 


J  MT-J; 


Minimum  J3  screen  coordinate 


IF,  I  ,  J 


three  elevation  points. 


J  MAX: 


Maximum  J3  screen  coordinate 


the  three  elevation  points 


r Ofi"F • 


The  frame  buffer. 


K ,  and  L 


Counters . 


I  FLAMES  : 


iotal  numoer 


'rHp'-  +- 


::ar,CL.a:  - 


constructed  fror 


-.E3  O  a  —  f- 


elevation  data  pci1 


VlN.  ivv.  .-.ivv  s.-. . 


.'/•.  ■/v  v  v.-’.1:/  v  v.v.v.v.  /.v., 


Copy  avivb-.jla  to  DTTCJ  does  ttct 

permit  fully  legible  reproduction 

13.  SUBROUTINE  FRAhE  FIL 

a.  Functions  performed 

Constructs  an  edge  list  from  the  three  transformed 
elevation  points.  These  edges  form  the  boundaries  f or  a 
polago^  fill  routine.  The  pixels  that  are  determined  to 
fail  u i t h i n  the  transformed  triangle  are  narked  in  the  frame 
buffer  . 


NODE  A,  "iCuE  B, 


and  NCDE  C 


j  m :  r-j ; 


Cutput 


“  PAVE  : 


-a.i.ng 


FILL  . 


The  numerical  designation  cf  the 


three  elevation  data  points  that 


make  a  triangular  panel  . 

Array  that  contains  the  IA  screen 


coordinates  of  the  transformed 


elevation  data  points. 

Array  that  contains  the  JA  screer 


coordinates  of  the  transformed 


elevation  data  points. 


Ma:<imum  JA  screen  coordinate  of  the 


three  elevation  data  points. 


M in : i mum  JA  screen  coordinate  of  the 


three  ele v alien  data  points. 


The  frame  buffer  arrsa  . 


y.v;.»  -•  -•>  -  ■  -• 


Called  routine 
None  . 


Routine  parameters 

I  I *JCFT :  T^s  cCw-'iiir'3ti3  zit'  tl 

maximum  point  of  a  lire. 


:clt3  : 


nn  ra  r  . 

OEL  : 

XNCDE1  and 
XNCDE  2: 
s  between  two 


The  amount  I  P  changes  for  a  o"o 
step  change  in  J p . 

The  total  JP  span  of  a  li~e. 
°rray  oontainmg  the  parameters 
the  three  lines  of  transformed 
triangular  panel  . 

Designates  the  number  of  ooc: 
Imes  for  a  gi' 


DX: 


jDE  : 


*  i  n  n  £ l t  and 
NCDE2 : 


Indicator  used  to  determine  the 
dirsotion  m  which  the  1 P  coor¬ 
dinates  are  counted. 

Pr  ray  containing  the  numerical 
designation  of  the  three  els'.at. 
data  ooints. 

Designators  for  deter tL-e 
e-evat;!”  point  that  oo-ta i "S  f 
highest  JQ  •  alue  between  t-o 
r o : "ts . 


AlV'jlw  AVViVyVxi  ■.yJ>  -V.V.  j  r_«V»iLd 


LHnrm«»»  •jwr’Trr.-yn  »■. »v»T»vr.'  yvT.'y. 

Co pf  available  lo  DTIC  goes  no! 
permit  fully  legibla  reproduction 


N  HIGH  and 


N  LOW: 


IT,  I U ,  IP 
ard  IS: 

J  B'JGKET: 


I i  and  I 2 : 


ICNT  and  LCNT : 


Used  with  UCEEI  and  UGGE2  Per 
determining  the  highest  JA  value 
between  t^c  pcirts. 

Counters , 

Array  that  contains  the  line  number 
designators  referencing  the  AEU  . 
Used  to  determine  the  number  of  I A 
coordinates  to  be  written  tc  the 
frame  buffer. 

Counters  fer  the  J  EUCKET . 
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APPENDIX  B 


3D-TRANSF0RMATI0N  PROGRAM  LISTING 
PROGRAM  TRAN_3_D 

THIS  PROGRAM  TAKES  AN  IMAGE  AND  ELEVATION  FILE  AND 
CONSTRUCTSTHE  REFERENCE  IMAGE  AND  ELEVATION  FILE. 
FROM  THESE  FILES  A  SYNTHESIZED  IMAGE  IS  PRODUCED. 

CHARACTER  ELFILE* 13 , IMFILE* 13 
BYTE  I MAGE (  512,512) 

INTEGER  IELEV2(  50,50), IA(  2500) , JA(  2500) 

INTEGER  I ROW, LROW, ICOL, LCOL, IENDN, IENDM, 

INTEGER  I  FRAME , JFRAME 

REAL  XI, Y1,Z1, FOCUS, DEPTH( 2500) 

REAL  A1 , A2 , B1 , B2 , Cl , C2 , OMEGA, PHI , KAPPA 

CALL  INPUT(  I ROW, LROW, ICCL, LCOL, ELFILE, IMFILE, 

I  FRAME, JFRAME) 

OPEN( UNIT=1 , FILE=ELFILE, STATUS= ' OLD ' ) 

OPEN(UNIT=2, FILE3 'ZFIL.DAT' , STATUS= ' NEW ' , 

ACCESS=' DIRECT' , RECORDSIZE=128 , MAXREC=5 12 ) 

OPEN( UN I T=3 , F I LE= ' XYZ . DAT' , STATUS=’NEW' , 

ACCESS3 ' SEQUENTIAL ’ , FORM3 1  FORMATTED '  ) 

OPEN( UNIT=4,FILE=IMFILE, STATUS3’ OLD' , ACCESS3 ' DIRECT ' , 
RECORDS I ZE=128,MAXREC=5 12) 

OPEN( UNIT=20, FILE3 'NODE.  DAT'  , STATUS3 'NEW'  , 

ACCESS3 'DIRECT' , RECORDSIZE=128, MAXREC=512 ) 

OPEN( UNIT=2 1 , FILE3 ' IMAGES. DAT' , STATUS3 'NEW' , 

ACCESS3 'DIRECT' , RECORDSI ZE=128 , MAXREC=5 12 ) 

CALL  TER_CROP( IROW, LROW, ICOL, LCOL, IELEV2 ) 

I ENDN=LROW- IROW+ 1 
I ENDM=LCOL- ICOL+ 1 
CALL  RE AD I MAGE (  IMAGE) 

CALL  TER_INTRP( IENDN, IENDM, IELEV2 ) 

CALL  REAL_EL(  IENDN, IENDM, IELEV2 ) 

CALL  TER_DMS(  IROW, LROW, ICOL, LCOL, IENDM) 

CALL  IM_REFI J(  I A, JA, I  FRAME , JFRAME , IENDM, IENDN) 

CALL  IM_REFAVG( I A, JA, IMAGE, IENDM, IENDN) 

CALL  AFFIN(  A1,A2,B1,B2,C1,C2) 

CALL  OBS_LOC(  XI, Y1,Z1, FOCUS) 

CALL  MEW_I J(  XI , Y1 , Z1 , FOCUS , A1 , A2 , B1 , B2 , Cl , C2 , 

IENDM, IENDN, IA,JA) 

CALL  NODE_DPTH(  XI, Y1,Z1, DEPTH, IENDM, IENDN) 

CALL  FILL(  IA, JA, DEPTH, IMAGE, IENDM, IENDN) 

CL0SE(  1) 

CL0SE(  2) 

CL0SE(  3 ) 

CLOSE( 4) 

CL0SE(  20) 
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CLOSE(  21) 


END 


SUBROUTINE  INPUT( IROW, LROW, ICOL, LCOL, ELFILE, IMFILE, 

I FRAME , JFRAME ) 


C 

C 

C 

c 

c 

c 

c 

c 

c 


I  ROW 
LROW 
ICOL 
LCOL 
ELFILE 
IMFILE 
I FRAME 
JFRAME 


INITIAL  ROW  OF  DESIRED  AREA 
LAST  ROW  OF  DESIRED  AREA 
INITIAL  COLUMN  OF  DESIRED  AREA 
LAST  COLUMN  OF  DESIRED  AREA 
ELEVATION  DATA  FILE  NAME 
IMAGE  DATA  FILE  NAME 
I  FRAME  NUMBER 
J  FRAME  NUMBER 
CHARACTER  ELFILE*13 , IMFILE* 13 
INTEGER  I ROW, LROW, ICOL, LCOL, IFRAME , JFRAME 


35 


45 


55 


WRITE(  6, *) 
WR I TE ( 6 , * ) 
READ( 5,35) 
WRITE( 6, * ) 
READ( 5,35) 
WRITE(  6, * ) 
READ( 5, 35) 
WRITE( 6, *) 
READ(  5,35) 
FORMAT ( 13) 
WRITE(  6,*) 
READ( 5,45) 
WRITE(  6,*) 
READ( 5, 45) 
FORMAT( A13 
WRITE( 6, *) 
READ( 5,55) 

write(  6,  * ) 

READ( 5,55) 
FORMAT ( 14) 
WRITE( 6, *) 
RETURN 
END 


' INPUT 
'  ENTER 


ELEVATION  AREA  DESIRED' 
MINIMUM  ROW  NUMBER  :  ' 


I  ROW 
' ENTER 
LROW 
'  ENTER 
ICOL 
'  ENTER 
LCOL 


MAXIMUM  ROW  NUMBER  :  ' 

MINIMUM  COLUMN  NUMBER  :  ' 

MAXIMUM  COLUMN  NUMBER  :  ' 


' INPUT 
ELFILE 


'  INPUT 
IMFILE 


THE  ELEVATION  DATA  FILE  NAME  :  ' 

THE  IMAGE  DATA  FILE  NAME  :  ' 


) 


'  INPUT 
IFRAME 


' INPUT 
JFRAME 


THE  I  FRAME  NUMBER  OF  IMAGE  :  ' 
THE  J  FRAME  NUMBER  OF  IMAGE  :  ' 


'*****  WAIT  APPROX.  1  MINUTE  FOR  SETUP*****' 


C 

C 


C 

C 

c 

c 

c 


★  ★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★TIT* 

SUBROUTINE  TER_CROP(  IROW, LROW, ICOL, LCOL, IELEV2 ) 


THIS  SUBROUTINE  READS  THE  ORIGINAL  TERRAIN  GRID 
AND  CONSTRUCTS  A  SMALLER  GRID  OF  THE  ELEVATION 
POINTS  DESIRED. 


CHARACTER  SWLAT*8, SWLON*8 , DELLAT*4, DELLON*4, 
CHARACTER  C0LS*4, R0WS*4 

INTEGER  IELEV1( 210,239) , IELEV2( 50,50) , IROW, LROW 
INTEGER  ICOL, LCOL, JV( 239) 

THE  VALUES  OF  JV, SWLON, SWLAT  DELLON, DELLAT, COLS 
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non 


AND  ROWS  ARE  IGNORED. 

READ(  1, 5)SWL0N, SWLAT, DELLON, DELLAT, COLS, ROWS 
FORMAT( 1X,2(A8,2X),4(A4,2X) ) 

DO  20  M=1 ,238 

READ(  1, 10) JV(M) ,  (  IELEV1( N, M) , N=1 , 210 ) 
FORMAT( 16, 2x, 201 6/( 8X, 2016 ) ) 

CONTINUE 

DO  40  N=IR0W, LROW 
DO  30  M=ICOL, LCOL 

IELEV2(  M+ 1- ICOL, N+ 1- IROW ) =IELEV1( N,  M ) 
CONTINUE 
CONTINUE 
RETURN 


SUBROUTINE  READ I MAGE ( IMAGE ) 

THIS  SUBROUTINE  READS  THE  IMAGE  DATA  INTO  AN  ARRAY. 
BYTE  IMAGE( 512,512) 

DO  10  IR=1 , 512 

READ( 4 , REC=IR ) ( I MAGE (  IC, IR) , IC=1,512) 

CONTINUE 

RETURN 

END 


SUBROUTINE  TER_INTRP(  IENDN, IENDM, IELEV2 ) 

THIS  SUBROUTINE  DETERMINES  THE  AVERAGE  VALUE  OF  THE 
ELEVATION  DATA  AND  ASSIGNS  THAT  VALUE  TO  ANY  UNKNOWN 
POINTS 

INTEGER  IENDN, IENDM, IELEV2( 50, 50) 

DATA  NEXT, I EL, IAVG/0,0,0/ 


DO  20  N=l, IENDN 
DO  10  M=l, IENDM 
IF(  IELEV2(  M,N).  EQ.  -32767)THEN 
GOTO  10 
ELSE 

NEXT=NEXT+ 1 
IEL=IEL+IELEV2(M,N) 

END  IF 
CONTINUE 
CONTINUE 
I AVG= I E  L/NEXT 

CHANGE  UNKNOWN  ELEVATION  VALUES  TO  THE 
CALCULATED  AVERAGE. 

DO  40  M=l, IENDN 
DO  30  M=l, IENDM 


IF( IELEV2( M,N) .  EQ.  -32767 )THEN 
IELEV2( M, N) =IAVG 
END  IF 

30  CONTINUE 
40  CONTINUE 
RETURN 
END 
C 

Q  **★★★★★******★★★★★★*★★★★★*★**★★★★★*★★**★★★*★*★★★★*★**•★*★★ 

SUBROUTINE  REAL_EL( IENDN, IENDM, IELEV2 ) 

C 

C  THIS  SUBROUTINE  CREATS  A  REAL  FILE  OF  THE  ELEVATIONS 

C 

INTEGER  IELEV2(  50,50), IENDN, IENDM 
REAL  AELEV( 239 ) 

C 

DO  20  N=l, IENDN 
K=0 

DO  10  M=1 , IENDM 
K=K+ 1 

AELEV( K)=IELEV2(M,N) 

10  CONTINUE 

WRITE( 2 , REC=N ) ( AELEV(  INT) , INT=1, IENDM) 

20  CONTINUE 

RETURN 
END 
C 

0  'k-k-k'k-k-k-k'k-k-k-k'k'k'k'k'k'k'k'k'k'k'k’k'k'k'k'k'k-k'k'k'k'k'k'k’k'k'k'k'k'k'k'k’k'k'k'k'k'k'k'k’kir'k'k'k'k'k 

SUBROUTINE  TER_DMS( IROW,LROW, ICOL,LCOL, IENDM) 

C 

C  THIS  SUBROUTINE  CONVERTS  EACH  ELEVATION  POINT  TO  ITS 

C  DMS  EQUIVALENT.  IT  USES  THE  FACT  THAT  EACH  ELEVATION 

C  POINT  REPRESENTS  A  ONE  SECOND  CHANGE  IN  LATITUDE 

C  OR  LONGITUDE  FROM  THE  NEXT  POINT. 

C  REFERENCE  POINT:  LAT. 037  DEG. : 22  MIN. : 47  SEC.  NORTH 

C  LON. -122  DEG. : 05  MIN. : 03  SEC.  WEST 

C  OUR  ELEVATION  POINTS  STAY  WITHIN  THE  BOUNDS  OF  037 

C  DEGREES  NORTH  AND  -122  DEGREES  WEST,  SO  THESE  VALUES 

C  WILL  BE  ASSIGNED. 

C 

IMPLICIT  DOUBLE  PRECISION  (A-Z) 

INTEGER  I ROW, LROW, I COL, LCOL, IENDM, IL, JL 
REAL  X , Y , 2 , LATD , LATM , L AT  S , AE  LEV (  2  3  9 ) 

REAL  LOND, LONM, LONS, HEIGHT 

PARAMETER( RLATM=22 .  ,RLATS=47.  ,RL0NM=5.  ,RL0NS=3.  ) 

C 

LATD=37. 

LOHD=- 122 . 

DO  -20  I  L=  I  ROW,  LROW 
K=IL/6Q 
LATM=RLATM+K 
LATS=RLATS+( IL-K*60) 

IF(  LATS.  GE.  60.  0 ) THEN 
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nnoooo  oo 


LATS=LATS-60. 0 
LATM=LATM+1. 0 
END  IF 

I l=IL-( IROW-1) 

READ( 2 , REC=I 1 ) ( AELEV(  INT)  , INT=1, IENDM) 

1=0 

DO  10  JL=IC0L/ LCOL 
K= JL/60 
L0NM=RL0NM-K 
L0NS=RL0NS-( JL-K*60) 

IF( LONS. LT. 0. 0 ) THEN 
L0NM=L0NM-1. 0 
LONS=LONS  +  60.  0 
END  IF 
LONM=-LONM 
LONS=-LONS 
1  =  1  +  1 

HEIGHT=AELEV( I ) 

CALL  DMS2XYZ(  LATD , LATM , LATS , LOND , LONM , LONS , 

HEIGHT, X, Y, Z ) 

WRITE(  3 , 99 )X, Y, Z 
99  FORMAT( 3( IX, FI 7. 7 )  ) 

10  CONTINUE 

20  CONTINUE 

ENDFILE( 3) 

REWIND( 3) 

RETURN 
END 

SUBROUTINE  IM_REFI J( IA, JA, I FRAME , JFRAME , IENDM, IENDN) 

THIS  SUBROUTINE  CONSTRUCTS  THE  I  AND  J  COORDINATE 
DATA  FOR  EACH  ELEVATION  POINT  AND  THEN  CONVERTS  THEM 
TO  SCREEN  COORDINATES  AND  STORES  THEM  IN  ARRAYS 
I A  AND  JA. 

IMPLICIT  DOUBLE  PRECISION  (A-Z) 

REAL  OMEGA, PHI , KAPPA, X, Y, Z , XI , Y1 , Z1 , XO , YO 
REAL  XIMA, YIMA  A1 , A2 , B1 , B2 , Cl , C2 , FOCUS 
INTEGER  I , J, IENDM, IENDN, IA( 2500) ,JA( 2500) 

INTEGER  I  FRAME, JFRAME, I TOT 

DATA  OMEGA, PHI , KAPPA/. 8341764, -. 4563699,3. 0761254/ 

DATA  X0,Y0/0. 000002,0. 0/ 

DATA  XI , Y1 , Zl/- 2 69 3 7 65. 9,-4304520. 4,3859018. 3/ 

DATA  Al, A2/20. 11323959,-6. 022849824/ 

DATA  El , B2/6. 016207940,20.  T0938801/ 

DATA  Cl , C2/-34954. 59484,-22566. 71593/ 

DATA  FOCUS/O. 153197/ 

I TOT= IENDN* IENDM 

do  io  :r=i,itot 
READ(  j,S,END=20)X,Y,Z 
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FORMAT( 3( IX, F17. 7)  ) 

CALL  PROJECT( X, Y, 2 , OMEGA , PHI , KAPPA, XO , YC , XI , Y1 , 21 , 
XIMA,YIMA, FOCUS) 

CALL  XY2 I J(  XIMA, YIMA, I , J, Al, A2 , Bl, B2 , Cl, C2 ) 

I A( I R ) = I - I FRAME 
JA( IR ) =4999- JFRAME- J 
10  CONTINUE 

20  REWIND(3) 

RETURN 
END 

★  ★*★★****★****★*★★★*•***★**★★★★*★*★★★★***★**★★★★★★*★*★**★*★ 
SUBROUTINE  DMS2XY2( LATD, LATM, LATS , LOND , LONM, LONS, 

LONS, HEIGHT, X, Y, 2) 

THIS  SUBROUTINE  CONVERTS  DMS  DATA  TO  X,Y,  AND  2 
GEORECTANGULAR  COORDINATES. 

IMPLICIT  DOUBLE  PRECISION  (A-2) 

REAL  PHI , LAMDA,N,X,Y, 2, LATD, LATM, LATS 
REAL  LOND, LONM, LONS, HEIGHT 
PARAMETER( PI=3.  14159265358793238) 

PARAMETER( Cl=180.  ,C2=60.  ,C3=3600.  ) 

PARAMETER(  E_SQUARE=0.  006768658 , A=6378206. 4) 

RADIAN=PI/C1 

PHI =(  LATD+LATM/C2  +  LATS/C3 ) *RADIAN 
LAMDA=( LOND+LONM/C2+LONS/C3 ) * RADIAN 
N=A/SQRT( 1-E_SQUARE*SIN(  PHI ) *SIN(  PHI ) ) 

X=( N+HE IGHT ) *COS( PHI )*COS( LAMDA ) 

Y=( N+ HEIGHT ) *COS(  PHI ) *SIN(  LAMDA) 

2=( N*( 1 -E_SQUARE ) +HE IGHT ) * S IN(  PHI ) 

RETURN 
END 

********************************************************* 

SUBROUTINE  PROJECT( X, Y, 2 , OMEGA, PHI , KAPPA, XO , YO , XI , 

Yl, 21, XIMA, YIMA, FOCUS) 

C 

C  THIS  SUBROUTINE  CONVERTS  THE  X, Y,2  GEORECTANGULAR 

C  COORDINATES  TO  XIMA  AND  YIMA  WHICH  ARE  IMAGE 

C  PLANE  COORDINATES. 

C 

IMPLICIT  DOUBLE  PRECISION  (A-2) 

REAL  Mil , M12 , M13 , M21 , M22 , M23 , M3  1 , M32 , M3 3 , DENOM 
REAL  X I MA , Y IMA , X , Y , 2 , OMEGA , PH I , KAPPA 
REAL  XO, Y0,X1, Yl, 21, FOCUS 
C 

Ml 1=C0S( PHI )*COS( KAPPA ) 

MI2=C0S(  OMEGA) *SIN(  KA?PA)  + 

5IN( OMEGA) *SIN( PHI)*COS( KAPPA) 

M13  =  SIi:(  OMEGA)  *SIN(  KAPPA )- 

CCS( OMEGA ) +SIN( PHI )*COS( KAPPA) 
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M21=-C0S( PHI )*SIN( KAPPA) 

M22=COS(  OMEGA) *COS(  KAPPA) - 

SIN( OMEGA) *SIN( PHI )*SIN( KAPPA) 

M23=SIN(  OMEGA ) *COS(  KAPPA)  + 

COS( OMEGA) *SIN( PHI )*SIN( KAPPA) 

M3 1=S IN( PHI ) 

M32=-SIN( OMEGA) *COS( PHI ) 

M33=COS( OMEGA) *COS( PHI ) 

DENOM=M3 1*(  X-Xl ) +M32*( Y-Yl ) +M33*( Z-Zl ) 

XIMA=XO-FOCUS*( Mll*( X-Xl )+M12*( Y-Yl ) +M13*(  Z-Zl ) ) /DENOM 

YIMA=YO-FOCUS*(  M21*(  X-Xl ) +M22* ( Y-Yl ) +  M23 *( Z-Zl ) ) /DENOM 

XIMA=XIMA* 1000000 

YIMA=Y IMA* 1000000 

RETURN 

END 


SUBROUTINE  XY2 I J( XIMA, YIMA, I , J, A1 , A2 , B1 , B2 , Cl , C2 ) 


THIS  SUBROUTINE  TAKES  THE  IMAGE  POINTS  XIMA, YIMA 
AND  CONVERTS  THEM  TO  I , J  ORIGINAL  IMAGE  COORDINATES. 


IMPLICIT  DOUBLE  PRECISION  (A-Z) 

REAL  XIMA, YIMA, A1 , A2 , Bi , B2 , Cl , C2 , DENOM 
INTEGER  I,J 


DENCM=A1*B2-BI*A2 

I=( ( XIMA- Cl ) *B2-( YIMA-C2)*B1) /DENOM 
J=-( ( XIMA- Cl ) *A2-( YIMA-C2 ) *A1 ) /DENOM 
RETURN 
END 


ir'k’k'k’k’ki'-k'k'k'k'k'k'k'kyr’k’k’k-k-k'k'k'k'k'k'k’k'k-k-k-kle'k-k’k'k'K'k-k-k'k-k'k'k’k'k  +  ’kyc'k'k'k'k'k’k'k 

SUBROUTINE  IM_REFAVG( IA, JA, IMAGE, IENDM, IENDN) 


THIS  SUBROUTINE  CONSTRUCTS  THE  FILE  THAT  IDENTIFIES 
THE  ELEVATION  POINTS  THAT  MAKE  UP  A  TRIANGULAR  PANEL 
AND  ITS  ASSOCIATED  SAMPLED  GREY  SCALE  VALUE. 


IMPLICIT  DOUBLE  PRECISION  (A-Z) 

REAL  SLOPE, YI NT 
BYTE  IMAGE( 512,512) 

INTEGER  I A ( 2500) , JA(  2500) , IY,M,N, I  GREY 1 , IGREY2 , L, LI 
INTEGER  IENDM, IENDN, NODE _ A , N0DE_B , N0DE_C , NCDE_D , IR 
INTEGER  I  COUNT 1 , I COUNTS , IT0T1, IT0T2 


1  90  I R=1 , I ENDN- 1 
IL=( IR-1)* IENDM 

1C  30  N=>IL,  :endk-i  +  il 

NODE  A=N 


:endm 


95 


. A  A  A  A  \ 


THREE  -  DIMENSIONAL  IHAQE  SERRATION  FROM  M  AERIAL 
PHOTOGRAPHS)  NAVAL  POSTGRADUATE  SCHOOL  HONTEREV  CA 
L  0  COLEMAN  SEP  07 


UNCLASSIFIED 


F/6  D/2 


SLOPE=(  JA(  NODE_C ) - JA( NODE_B ) ) * 1 . 0/ 

(  IA(NODE_C)-IA(NODE_B)  )*1.  0 
YINT=1.  0* JA(  NODE_B ) -SLOPE* IA( NODE_B ) 

I TOT 1=0 
IT0T2=0 
ICOUNT1=0 
ICOUMT2=0 

DO  70  M=IA(  NODE_B ) , IA( NODE„A) 

IY=(  SLOPE*M+YINT) 

DO  50  L=JA(NODE_B) / IY 
IT0T1=IT0T1+ IMAGE( M, L) 

IC0UNT1=IC0UNT1+1 

CONTINUE 

IF(M.  LT.  IA(NODE_A)  )THEN 
DO  60  L=IY+ 1 , JA( NODE_C ) 

ITOT2=ITOT2+IMAGE( M,  L) 

I COUNT2 = I COUNT2 + 1 
CONTINUE 
END  IF 
CONTINUE 
Ll=( N-( IR-1) )*2 
IGREY1=IT0T1/IC0UNT1 
I GRE Y2 = I TOT2 / 1 COUNT2 
NODE_D=NODE_C+l 

WRITE(  20,REC=L1-1)NODE_A,NODE_B,NODE_C, I GREY 1 
WRITE(  20,REC=L1)NODE_B/NODE_D,NODE_C, I GREY 2 
CONTINUE 
CONTINUE 
RETURN 
END 

**************************************************** 

SUBROUTINE  EXT_ORIEN(  Mil , M12 , M13 , M2 1 , M22 , M23 , 

M31,M32,M33, I ENDM , I ENDN ) 

THIS  ROUTINE  DETERMINES  THE  M  MATRIX  PARAMETERS 
FOR  ROTATION  OF  THE  IMAGE  PLANE  TO  DESIRED 
LOCATION  FOR  VIEWING  IN  THE  SYNTHESIZED  IMAGE. 

IMPLICIT  DOUBLE  PRECISION  (A-Z) 

REAL  MAGN_X , MAGN_Y , MAGN_Z , X , Y , Z 

REAL  X_C0RD(  2500 ) , Y_CORD( 2500 ) , Z_C0RD( 2500) 

REAL  X_VECX , X_VECY , X_VECZ , Y_VECX, Y_VECY, Y_VECZ 
REAL  Z_VECX, Z  VECY,Z_VECZ 
INTEGER  IENDM, IENDN, ITOT, IP, IR 

I TOT= I ENDM* I ENDN 
DO  10  IR=1 , ITOT 
READ( 3,5, END=20 )X, Y, Z 
X_C0RD( IR)=X 
Y_C0RD( IR)=Y 
Z_CORD( I R ) =Z 
FORMAT {  3(  IX, F17.  7) ) 


10  CONTINUE 
20  REWIND( 3 ) 

Y_VECX=-(  X_C0RD( 1 ) ) 

Y_VECY=-(  Y_CORD( 1 ) ) 

Y_VECZ=-(  Z_CORD( 1 ) ) 

IP=ITOT- IENDM+1 
Z_VECX=X_CORD( 1 ) -X_CORD( IP) 

Z_VECY=Y_CORD( l)-Y_CORD(  IP) 

Z_VECZ=Z_CORD( 1 ) -Z_CORD( IP ) 

C  USE  THE  CROSS  PRODUCT  OF  Y  CROSS  Z  TO  OBTAIN 

C  THE  X  VECTOR. 

X_VECX=(  (  Y_VECY*Z_VECZ)-( Y_VECZ*Z_VECY) ) 

X_VECY=(  (  Y_VECZ*Z_VECX) -(  Y_VECX*Z_VECZ ) ) 

X_VECZ=(  (  Y_VECX*Z_VECY)-( Y_VECY*Z_VECX) ) 

MAGN_Z=SQRT(  (  Z„VECX**2 ) +(  Z_VECY**2)  +  (  Z_VECZ**2) ) 
MAGN_X=SQRT{  (  X_VECX**2 )  +  ( X_VECY**2 )  +  ( X_VECZ**2 ) ) 
MAGN_Y=SQRT(  (  Y_VECX**2)  +  (  Y_VECY**2)+( Y_VECZ**2) ) 

Ml 1=X_VECX/MAGN_X 
M 1 2  =X_VECY/MAGN_X 
M 1 3  =X_VECZ/MAGN_X 
M2 1=Y_VECX/MAGN_Y 
M2  2 =Y_VECY/MAGN_Y 
M2 3=Y_VECZ/MAGN_Y 
M3 1=Z_VECX/MAGN_Z 
M3 2 =Z_VEC Y/MAGN_Z 
M33=Z_VECZ/MAGN_Z 
RETURN 
END 
C 

C  A********************************************************* 

SUBROUTINE  AFFIN( A1 , A2 , B1 , B2 , Cl , C2 ) 

C 

C  THIS  SUBROUTINE  ASSIGNS  OR  CALCULATES  THE 

C  COEFFICIENTS  TO  BE  UTILIZED  IN  THE  TRANSFORM  FROM 

C  IMAGE  COORDINATES  TO  SCREEN  COORDINATES. 

C 

IMPLICIT  DOUBLE  PRECIS ION( A-Z) 

REAL  Al , A2 , B1 , B2 , Cl , C2 , XIMA_MAX, YIMA_MAX, I_MAX, J_MAX 
DATA  I_MAX, J_MAX/512.  0,512.  0/ 

C 

XIMA_MAX=1600.  0 
YIMA_MAX=1600.  0 
C1=XIMA_MAX 
C2=0.  0 
A2=0. 0 
B1=0. 0 

Al=-XIMA_MAX/( I_MAX* 1. 0) 

B2=YIMA_MAX/(  J_MAX*1.  0) 

RETURN 

END 

C 

0  ********************************************************** 
SUBROUTINE  OBS_LOC{  XI , Y1 , Z1 , FOCUS ) 
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c 

C  THIS  SUBROUTINE  CALCULATES  THE  NEW  OBSERVER  X1,YL,Z1 

C  LOCATION  FROM  DESIRED  LAT.  AND  LONG.  INPUTS  AS  WELL 

C  AS  PROVIDE  THE  FOCAL  LENGTH. 

C 

IMPLICIT  DOUBLE  PRECISION  (A-Z) 

REAL  LATD , LATM , LATS , LOND , LONM , LONS , HE I GHT 
REAL  X1,Y1,Z1,X,Y,Z, FOCUS 
C 

LATD=37.  0 

WRITE( 6,*) ' INPUT  OBSERVER  LATITUDE  IN-MINUTES( REAL) : ' 
READ( 5,5) LATM 

WRITE( 6 , * ) '  -SECONDS( REAL): ' 

READ( 5, 5) LATS 
LOND=-122.  0 

WRITE( 6,*) ' INPUT  OBSERVER  LONGITUDE  IN-MINUTES( REAL) : ' 
READ( 5,5) LONM 

WRITE( 6 ,  *  )  '  -SECONDS (REAL): ' 

READ( 5,5) LONS 
5  FORMAT( F5. 1 ) 

WRITE( 6,*) ' INPUT  OBSERVER  HEIGHT-METERS(  REAL) :  ' 

READ( 5, 10) HEIGHT 
10  F0RMAT(F6.  1) 

FOCUS=0. 015 

CALL  DMS2XYZ(  LATD, LATM, LATS, LOND, LONM, LONS, HEIGHT, 

X, Y, Z ) 

X1=X 


Y1=Y 

Z1=Z 

RETURN 

END 


C 


c 

c 

c 

c 


10 


****************************************************** 
SUBROUTINE  NEW_I J( XI , Y1 , Z1 , FOCUS , A1 , A2 , B1 , B2 , Cl , C2 , 

IENDM, IENDN, IA, JA) 


THIS  SUBROUTINE  COMPUTES  THE  NEW  I A  AND  JA  SCREEN 
COORDINATES  FROM  THE  GIVEN  OBSERVER  LOCATION. 


IMPLICIT  DOUBLE  PRECISION  (A-Z) 

REAL  XI , Y1 , Z1 , FOCUS , Mil , M12 , M13 , M21 , M22 , M23 , M31 
REAL  XO, YO , X, Y, Z, XIMA, YIMA, A1 , A2 , B1 , B2 , Cl , C2 
REAL  M3  2 , M3  3 , DENOM 

INTEGER  IENDM, IENDN, ITOT, IR, IA( 2500) ,JA( 2500) , I, J 
DATA  X0,Y0/0. 0,0.  0/ 


I TOT= IENDN* IENDM 
DO  10  IR=1 , 2500 
IA( IR)=0 
JA( IR) =0 
CONTINUE 

CALL  EXT_0RIEN( Mil , M12 , M13 , M21 , M22 , M23 , M31 , M32 , M33 , 
IENDM, IENDN) 
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DO  20  IR=1 , ITOT 
READ( 3,15, END=30 )X, Y, Z 
15  FORMAT( 3( IX, F17. 7 ) ) 

DEN0M=M31*(X-X1)+M32*( Y-Y1)+M33*( Z-Zl) 

XIMA=XO- FOCUS *(  Mll*(  X-Xl ) +M12*( Y-Yl )+M13*( Z-Zl ) )/ 
DENOM 

YIMA=YO-FOCUS*(  M21*{ X-Xl  )+M22*(  Y-Yl )+M23*(  Z-Zl)  )/ 
DENOM 

XIMA=XIMA* 1000000.  0 
YIMA=Y IMA* 1000000. 0 

CALL  XY2 1 J(  XIMA, YIMA, I , J, Al, A2 , B1 , B2 , Cl, C2 ) 

IA( IR)=I 
JA( IR)=J 
20  CONTINUE 

30  REWIND( 3 ) 

RETURN 

END 

C 

SUBROUTINE  NODE_DPTH( XI , Y1 , Z1 , DEPTH, IENDM, IENDN) 

C 

C  THIS  SUBROUTINE  CALCULATES  THE  DISTANCE  OF  THE 

C  TRANSLATED  ELEVATION  POINTS  TO  THE  NEW  OBSERVER 

C  LOCATION. 

C 

IMPLICIT  DOUBLE  PRECI SION( A-Z ) 

REAL  X1,Y1,Z1,DEPTH( 2500) 

INTEGER  IENDM, IENDN, IR, ITOT 
C 

I TOT= IENDM* IENDN 
DO  10  IR=l,ITOT 
READ( 3,5, END=20 )X, Y, Z 
5  FORMAT( 3( IX, F17.  7 ) ) 

DEPTH( IR)=SQRT(  (  (  Xl-X)**2 )  +  (  (  Yl-Y)**2 )  +  (  (  Zl-Z)**2) ) 

10  CONTINUE 

20  REWIND( 3 ) 

RETURN 

END 

C 

SUBROUTINE  FILL( IA, JA, DEPTH, IMAGE, IENDM, IENDN) 

C 

C  THIS  SUBROUTINE  DETERMINES  THE  HIDDEN  SURFACES  AND 

C  CONSTRUCTS  THE  TRANSLATED  IMAGE.  A  Z_BUFFER  IS  USED 

C  TO  HOLD  THE  COMPUTED  DEPTHS  FROM  THE  OBSERVER  TO  THE 

C  GEOGRAPHIC  POSITION.  A  PLANE  EQUATION  IS  CONSTRUCTED 

C  FROM  THREE  ELEV.  POINTS.  THIS  EQUATION  IS  USED  TO 

C  DETERMINE  THE  DEPTH  OF  ALL  POINTS  WITHIN  THE  PLANE. 

C 

IMPLICIT  DOUBLE  PRECISION A-Z) 

BYTE  IMAGE( 512,512) 

REAL  DEPTH(  2500 ) , Z_BUFF(  5 12 , 512 ) , XI , X2 , X3 , Y1 , Y2 , Y3 
REAL  Z1 , Z2 , Z3 , Z_DPTH, Cll , C12 , C13 , C21 , C22 , C23 , C31 
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REAL  C3  2 , C3  3 , DET , A_COEF , B_COEF , C.COEF 
INTEGER  I GREY , NODE_A , NODE_B , NODE_C , I_MIN, I_MAX, J_MIN 
INTEGER  J_MAX, FRAME ( 512,512 ) , IA( 2500) , JA( 2500) , IENDM 
INTEGER  IENDN, IR,  I ,  J,  K,  L, I PLANES 

FIRST  DETERMINE  THE  NUMBER  OF  PLANES  AND  INITIALIZE 
THE  Z_BUFF£R  AND  IMAGE  TO  0. 

IPLANES=(  (  IENDM-1 ) *2 ) *(  IENDN-1) 

DO  20  K=l, 512 
DO  10  L=1 , 512 
I MAGE ( L, K)=0 
Z_BUFF( L, K) =0 
CONTINUE 
CONTINUE 

DETERMINE  THE  COEFFICIENTS  OF  THE  PLANE  EQUATION 
FROM  THE  THREE  ELEVATION  POINTS. 

DO  70  IR=1 , IPLANES 

READ( 20 , REC=IR) NODE_A , NODE_B , NODE_C , I GREY 
X1=IA( NODE_A ) 

Y1=JA( NODE_A) 

Z1=DEPTH( NODE_A ) 

X2=IA(NODE_B) 

Y2=JA( NODE_B) 

Z2=DEPTH( NODE_B) 

X3=IA( NODE_C) 

Y3=JA( NODE_C) 

Z3=DEPTH( NODE_C) 

DETERMINE  THE  COFACTOR  ELEMENTS 
Cll=(  (  Y2*Z3)-( Y3*Z2) ) 

C12=-( (X2*Z3)-(X3*Z2)) 

C13=(  (X2*Y3)-(X3*Y2)  ) 

C2 1=-(  (Y1*Z3)-(Y3*Z1)  ) 

C22=(  (X1*Z3)-(X3*Z1)  ) 

C23=-(  (X1*Y3)-(X3*Y1)  ) 

C31=(  ( Yl*Z2)-( Y2*Z1) ) 

C32=-(  (X1*Z2)-(X2*Z1)  ) 

C33=(  (  X1*Y2 ) -( X2*Y1 )  ) 

CALCULATE  THE  DETERMINANT 

DET=( X1*C11 )+( Y1*C12 ) +( Z1*C13 ) 

THE  COEFFICIENTS  ARE  DETERMINED  FROM  MULTIPLYING 
THE  reciprocal  of  the  determenant  with  the 
transpose  of  the  cofactors  called  the  adjoint. 

A__COEF=-(  ( C11+C21+C31  )/DET) 

B_COEF=-(  ( C12+C22+C32 )/DET) 

C_COEF=-( ( C13+C23+C33 ) /DET ) 

IF(  C_C0EF.  EQ.  0.  0)GOTO  70 

DETERMINE  THE  MAXIMUM  AND  MINIMUM  VALUES  TO  TEST 


C  FOR  THE  FILL  ALGORITHUM. 

C 

I_MAX=MAX( IA( NODE_A ) , IA( NODE_B) , IA( NODE_C) ) 

I_MIN=MIN(  IA( NODE_A) , IA(  NODE_B ) , IA(  NODE_C) ) 

J_MAX=MAX(  JA(  NODE_A ) , JA( NODE_B) , JA( NODE_C) ) 

J_MIN=MIN( JA(  NODE_A ) , JA(  NODE_B ) , JA(  NODE_C ) ) 

IF(  I_MIN.  LT.  1.  OR.  I_MAX.  GT.  512)GOTO  70 
IF(  J_MIN.  LT.  1.  OR.  J_MAX.  GT.  512  )GOTO  70 
C 

C  CLEAR  THE  REFERENCE  FRAME  BUFFER  AND  CALL  THE 

C  FRAME  FILL  SUBROUTINE. 

C 

DO  40  L=I_MIN, I_MAX 
DO  30  K=J_MIN, J_MAX 
FRAME ( L/K)=0 
30  CONTINUE 

40  CONTINUE 

CALL  FRAME_FIL(NODE_A,NODE_B,NODE_C, FRAME, I A, JA, 

J_MAX , J_MIN ) 

DO  60  J=J_MIN, J_MAX 
DO  50  I=I_MIN, I_MAX 
IF( FRAME (  I, J).  EQ.  1 ) THEN 

Z_DPTH=-(  l  +  (  A_COEF*I )  +  (  B_COEF*J) )/C_COEF 
IF(  Z_BUFF(  I ,  J).  EQ.  0.  0.  OR.  Z_DPTH.  LT.  Z_BUFF(  I,  J)  )THEN 
Z_BUFF( I, J)=Z_DPTH 
I MAGE ( I , J ) = I GREY 

end  if 

END  IF 

50  CONTINUE 

60  CONTINUE 

70  CONTINUE 

DO  90  J=l, 512 

WRITE( 2 1 , REC=J ) ( I MAGE ( I , J) , 1=1 , 512 ) 

90  CONTINUE 

RETURN 
END 
C 


SUBROUTINE  FRAME_FIL(  NODE_A , NODE_B , NODE_C , FRAME, 

IA, JA, J_MAX, J_M I N ) 

THIS  SUBROUTINE  CONSTRUCTS  AN  EDGE  LIST  FROM  THE 
THREE  NODES  PASSED  IN  THE  ROUTINE.  THESE  EDGES 
ARE  USED  IN  A  POLYGON  FILL  ROUTINE  USEING  A  FRAME 
BUFFER  AND  Y.BUCKET 
IMPLICIT  DOUBLE  PRECI S ION( A-Z ) 

REAL  I_INCPT, DELTA_I , DELTA_J, AEL(  3,4) 

REAL  XNCDE1 , XN0DE2 , DX 

INTEGER  NODE  A, NODE  B , NODE_C, FRAME( 5 12 , 512 ) , I A( 2500 ) 
INTEGER  JA(  2500 ) , J_MAX, J_MIN, NODE(  4 ) , IS , N0DE1 , N0DE2 
INTEGER  N_H I GH , N_L0W , IT, IU, J_BUCKET(  512) , IR, II, 12 
INTEGER  ICNT, LCNT 


1 


a 


4 

i 


a 


DO  10  IS=1 , 512 
J_BUCKET(  IS)=0 
CONTINUE 
DO  30  IS=1 , 3 
DO  20  IR=1 , 4 
AEL(  IS,  IR)  =0.  0 
CONTINUE 
CONTINUE 
NODE( l)=NODE_A 
NODE( 2 ) =NODE_B 
NODE( 3 )=NODE_C 
NODE( 4 ) =NODE_A 
DO  40  IS=1 , 3 
NODEl=NODE( IS) 

NODE2=NODE(  IS+1) 

IF( JA( NODE1). GE. JA( NODE2 ) ) THEN 
N_H I GH=NODE 1 
N_LOW=NODE2 
ELSE 

N_HIGH=NODE2 
N_LOW=NODE 1 
END  IF 

I_INCPT=IA(N_HIGH) 

DELTA_J=( JA( N_HIGH) - JA( N_LOW) ) 

IF(DELTA_J.  EQ.  0.  0) THEN 
GOTO  40 
ELSE 

DELTA_I=-( IA( N_HIGH)-IA( N_LOW) )/DELTA_J 
END  IF 

IT=JA( N_H I GH ) 

IU=J_BUCKET(  IT) 

IF(  IU.  EQ.  0  ) THEN 
J_BUCKET(  IT)  =  IS 
ELSE 

AEL( IU, 4)=IS 
END  IF 

AEL( IS, 1 )=I_INCPT 
AEL( IS, 2 )=DELTA_I 
AEL( IS, 3 )=DELTA_J 
CONTINUE 

IT=J_BUCKET(  J_MAX) 

IU=INT( AEL(  IT, 4) ) 

XN0DE1=AEL(  IT, 1) 

XNODE2=AEL(  IU, 1) 

DX=XN0DE1-XN0DE2 
IF(  DX.  LE.  0.  0 ) THEN 
I 1=NINT( XNODE1 ) 

I2=NINT( XNODE2 ) 

ELSE 

I 1=NIMT( XNODE2 ) 

I2=NINT( XN0DE1 ) 

END  IF 

DO  50  IR=I 1 , 12 


FRAME ( IR, J_MAX) =1 
CONTINUE 

AEL(  IT,  3  )  =AEL(  IT,3)-1.  0 
AEL( IU, 3 )=AEL(  IU,3)-1.  0 
I CNT= J_MAX- 1 
LCNT=J_MIN 

DO  70  IS=ICNT, LCNT, -1 
IF( J_BUCKET( I S ) . EQ. 0 ) THEN 
XN0DE1=XN0DE1+AEL( IT, 2 ) 
XNODE2=XNODE2+AEL( IU, 2 ) 

ELSE  IF(  AEL(  IT,  3  ) .  LE.  0.  0 ) THEN 
IT=J_BUCKET( IS) 
XN0DE1=AEL( IT, 1) 

XNODE2  =XNODE2  + AEL( IU, 2 ) 
ELSE 

IU=J_BUCKET( IS) 
XN0DE1=XN0DE1+AEL( IT, 2) 
XNODE2=AEL( IU, 1) 

END  IF 

DX=XNODEl-XNODE2 
I F(  DX.  LE.  0.  0 ) THEN 
I 1=NINT( XN0DE1 ) 

I2=NINT( XNODE2 ) 

ELSE 

I 1=NINT( XNODE2 ) 

I 2=NINT( XNODE1 ) 

END  IF 

DO  60  IR=I1 , 12 
FRAME ( IR, IS)=1 
CONTINUE 

AEL(  IT,  3  )  =AEL{  IT,3)-1.  0 
AEL(  IU,  3  )=AEL(  IU,3)-1.  0 
CONTINUE 
RETURN 
END 
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